/* Author : 
 *    Daohang Shi 
 *    dshi7@wisc.edu 
 *    Fall  2013  
 *    */

#include "CGR.h"

#define DELTA_EDGE_WEIGHT 0.000001

bool  compare_pair_int_double ( const pair<int, double> &p1, const pair<int, double> &p2 ) { return  p1.second > p2.second; }

bool  compare_pair_int_GridBox ( const pair<int, GridBox > &p1, const pair<int, GridBox> &p2 ) { return  p1.second.dx * p1.second.dy > p2.second.dx * p2.second.dy; }

bool  compare_pair_GridBox_double ( const pair<GridBox, double> &p1, const pair<GridBox, double> &p2 ) { return p1.second>p2.second; }

bool  compare_pair_string_double ( const pair<string, double> &p1, const pair<string, double> &p2 ) { return  p1.second < p2.second; }

bool  compare_netroute_with_smaller_weight ( const pair<NetRoute, double> &p1, const pair<NetRoute, double> &p2 )  { return  p1.second < p2.second; }

bool  compare_subnet_with_smaller_weight ( const pair<TwoTerCoarseNet, double> &p1, const pair<TwoTerCoarseNet, double> &p2 )  { return  p1.second < p2.second; }

struct more_of : public binary_function<pair<uint32_t,int> , pair<uint32_t,int> , bool> {
  bool operator()(const pair<uint32_t,int> a, const pair<uint32_t,int> b){
    return ( a.second < b.second);
  }
};

int   CGR::xyzToGridNetTerminalId ( int x, int y, int z ) 
{
  if ( x<0 || x >= fullGridX )
    return  -1;
  if ( y<0 || y >= fullGridY )
    return  -1;
  if ( z<0 || z >= numLayers )
    return  -1;
  return  ( z*fullGridY*fullGridX + y*fullGridX + x );
}

/*! \fn CoarseNet CGR::getModifiedNetInBox ( CoarseNet net )
  \brief Read net and return a modified net inside box region
  \param CoarseNet net
  \returns CoarseNet
  */
CoarseNet CGR::genSingleModifiedNetInBox ( CoarseNet &net, GridBox &box ) 
{
  /* create a new net data structure */
  CoarseNet new_net;
  CoarseNet &modified_net = new_net;
  modified_net.old_index = net.old_index;
  /* update vector<Segment> net.sg and set<int> net.terminal_pins */ 
  for ( size_t  i=0;  i<net.sg.size(); i++ ) {
    Segment &seg = net.sg[i];
    addModifiedSegmentInNet ( seg, modified_net, box );
  }
  //  if ( modified_net.terminal_pins.size() > 0 )
  if ( modified_net.inside_terminal_pins.size() > 0 || modified_net.boundary_terminal_pins.size() > 0 ) {
    /* merge inside pins and boundary pins */
    modified_net.terminal_pins.clear();
    modified_net.terminal_pins.insert ( modified_net.inside_terminal_pins.begin(), modified_net.inside_terminal_pins.end() );
    modified_net.terminal_pins.insert ( modified_net.boundary_terminal_pins.begin(), modified_net.boundary_terminal_pins.end() );
  }
  return  new_net;
}

/*! \fn   void    CGR::addModifiedSegmentInNet ( Segment seg, CoarseNet &bnet ) 
  \brief Read segment and insert all the inbox segments and edges into the net object
 ************  ALSO *********** add vertex into net
 \param Segment &seg CoarseNet &bnet
 */
void    CGR::addModifiedSegmentInNet ( Segment &seg, CoarseNet &modified_net, GridBox &box ) 
{
  /* read each segment in old net class
   *  determine if part of it lays inside bounding box 
   *    by checking if two ending points lay inside box */
  int s = seg.s;
  int t = seg.t;
  int ending_point_on_boundary;
  CoarseNet &old_net = nets[map_from_net_old_index_to_idnex_in_nets[modified_net.old_index]];

  /* this segment may probably pass through the box region */
  if ( !isInBox (s,box) && !isInBox (t,box) ) {
    if ( seg.type=="VIA" )
      return;
    if ( seg.type=="HORIZ" ) {
      if ( V[s].x < box.x0 && V[t].x >= box.x0+box.dx && V[s].y>=box.y0 && V[s].y<box.y0+box.dy ) {
        assert ( V[s].y==V[t].y );
        int s1 = xyzToGridNetTerminalId ( box.x0, V[s].y, V[s].z );
        int t1 = xyzToGridNetTerminalId ( box.x0+box.dx-1, V[s].y, V[s].z );
        Segment bseg;
        bseg.s = s1;  bseg.t = t1;  bseg.type = "HORIZ";
        modified_net.sg.push_back ( bseg );
        modified_net.boundary_terminal_pins.insert ( s1 );
        modified_net.boundary_terminal_pins.insert ( t1 );
      }
      else
        return;
    }
    if ( seg.type=="VERT" ) {
      if ( V[s].y < box.y0 && V[t].y >= box.y0+box.dy && V[s].x>=box.x0 && V[s].x<box.x0+box.dx ) {
        assert ( V[s].x==V[t].x );
        int s1 = xyzToGridNetTerminalId ( V[s].x, box.y0, V[s].z );
        int t1 = xyzToGridNetTerminalId ( V[s].x, box.y0+box.dy-1, V[s].z );
        Segment bseg;
        bseg.s = s1;  bseg.t = t1;  bseg.type = "VERT";
        modified_net.sg.push_back ( bseg );
        modified_net.boundary_terminal_pins.insert ( s1 );
        modified_net.boundary_terminal_pins.insert ( t1 );
      }
      else
        return;
    }
  }

  if ( isInBox (s,box) && isInBox (t,box) ) {
    vector<int> node_vector = getNodeVectorFromSegment ( seg );
    for ( vector<int>::iterator itr_node = node_vector.begin(); itr_node != node_vector.end(); itr_node ++ ){
      if ( isOnBoundary ( *itr_node, box ) )
        modified_net.boundary_terminal_pins.insert ( *itr_node );
      else {
        if ( old_net.terminal_pins.find(*itr_node) != old_net.terminal_pins.end() )
          modified_net.inside_terminal_pins.insert ( *itr_node );
      }
    }
    modified_net.sg.push_back ( seg );
  }

  if ( isInBox (s,box) && !isInBox (t,box) ) {
    if ( seg.type == "HORIZ" ) 
      ending_point_on_boundary = xyzToGridNetTerminalId ( box.x0+box.dx-1, V[s].y, V[s].z );
    if ( seg.type == "VERT") 
      ending_point_on_boundary = xyzToGridNetTerminalId ( V[s].x, box.y0+box.dy-1, V[s].z );

    if ( ending_point_on_boundary== s )
      return;

    modified_net.boundary_terminal_pins.insert ( ending_point_on_boundary );
    //  s
    if ( isOnBoundary ( s, box ) )
      modified_net.boundary_terminal_pins.insert ( s );
    else {
      if ( old_net.terminal_pins.find(s) != old_net.terminal_pins.end() ) 
        modified_net.inside_terminal_pins.insert ( s );
    }
    //  bseg
    Segment bseg;
    bseg.s = s;  bseg.t = ending_point_on_boundary;  bseg.type = seg.type;
    modified_net.sg.push_back ( bseg );
  }
  if ( !isInBox (s,box) && isInBox (t,box) ) {
    if ( seg.type == "HORIZ" ) 
      ending_point_on_boundary = xyzToGridNetTerminalId ( box.x0, V[t].y, V[t].z );
    if ( seg.type == "VERT" ) 
      ending_point_on_boundary = xyzToGridNetTerminalId ( V[t].x, box.y0, V[t].z );

    if ( ending_point_on_boundary== t )
      return;

    modified_net.boundary_terminal_pins.insert ( ending_point_on_boundary );
    //  t
    if ( isOnBoundary ( t, box ) )
      modified_net.boundary_terminal_pins.insert ( t );
    else {
      if ( old_net.terminal_pins.find(t) != old_net.terminal_pins.end() ) 
        modified_net.inside_terminal_pins.insert ( t );
    }
    //  bseg
    Segment bseg;
    bseg.s = ending_point_on_boundary;  bseg.t = t;  bseg.type = seg.type;
    modified_net.sg.push_back ( bseg );
  }
}

/*! \fn   void    CGR::insertEdgesFromSegments ( Segment &seg, CoarseNet &bnet ) 
  \brief insert all the edges in segment into the given net
  \param Segment &seg CoarseNet &bnet
  */
void    CGR::insertEdgesFromSegments ( Segment &seg, CoarseNet &bnet ) 
{
  //  InsertEdgesFromSegmentsToNet
  vector<int> edge_vector = getEdgeVectorFromSegment ( seg );
  for ( size_t  i=0; i< edge_vector.size(); i++ )
    bnet.edges.insert ( edge_vector[i] );
}

/* 
 * Net decomposition determines both the pin ordering and the edges in each subnet 
 * */
void  CGR::netDecomposition ( CoarseNet &net ) 
{
  //  Using values : 
  //    set<int> terminal_pins 
  //    set<int> edges
  //  Assigned values :
  //    vector<TwoTerCoarseNet> subnets
  /* construct a graph in adjacent form during reading the segment vector */
  map<int, set<int> > tree_in_adjacent_form;  //  the value for each key is the all its adjacent nodes
  for ( set<int>::iterator itr_edge = net.edges.begin();itr_edge != net.edges.end(); itr_edge++ ) {
    int s = E[*itr_edge].srcPoint;
    int t = E[*itr_edge].snkPoint;
    if ( tree_in_adjacent_form.count (s) != 0 )  {  //  already exists 
      tree_in_adjacent_form[s].insert ( t );
    }
    else {
      set<int>  myset;
      myset.insert ( t );
      tree_in_adjacent_form.insert ( make_pair ( s, myset ) );
    }
    if ( tree_in_adjacent_form.count (t) != 0 )  {  //  already exists 
      tree_in_adjacent_form[t].insert ( s );
    }
    else {
      set<int>  myset;
      myset.insert ( s );
      tree_in_adjacent_form.insert ( make_pair ( t, myset ) );
    }
  }

  /* pin ordering */
  vector<int> terminal_pins ( net.terminal_pins.begin(), net.terminal_pins.end() );
  set<int>    visited_terminal_pins;

  //  number of subnets is terminal_pins.size()-1
  net.subnets.clear();
  int current_terminal_pin = terminal_pins[0];
  int next_terminal_pin;
  for ( int num_iteration = 0; num_iteration < terminal_pins.size()-1;  num_iteration++ ) {
    visited_terminal_pins.insert ( current_terminal_pin );
    set<int>  path;
    set<int>  &ref_path = path;
    next_terminal_pin = getNearestUnvisitedTerminalPin ( terminal_pins, tree_in_adjacent_form, current_terminal_pin, visited_terminal_pins, ref_path ); 
    //  define a new subnet
    TwoTerCoarseNet subnet;
    subnet.s = current_terminal_pin;
    subnet.t = next_terminal_pin;
    subnet.edges = path;
    //  
    current_terminal_pin = next_terminal_pin;
    //  add a new subnet
    net.subnets.push_back ( subnet );
  }
  //  order all the subnets according to their subnet weights
  vector<pair<TwoTerCoarseNet, double> >  subnet_weights_pair;
  for ( vector<TwoTerCoarseNet>::iterator it = net.subnets.begin(); it != net.subnets.end(); it++ )
    subnet_weights_pair.push_back ( make_pair(*it, computeSubnetEdgeWeight(*it)) );
  sort ( subnet_weights_pair.begin(), subnet_weights_pair.end(), compare_subnet_with_smaller_weight );
  net.subnets.clear();
  for ( vector<pair<TwoTerCoarseNet,double> >::iterator it = subnet_weights_pair.begin(); it != subnet_weights_pair.end(); it++ )
    net.subnets.push_back ( (*it).first );
}

/* find congested edges among vector<int> E_2d */
void  CGR::findCongestedEdges () 
{
  congested_edges.clear();
  for ( vector<int>::iterator it = E_id.begin(); it != E_id.end(); it++ ) 
  {
    if ( E[*it].cap!=0 && E[*it].type!="VIA"  && E[*it].congestion_ratio > TCR ) {         
      congested_edges.push_back ( make_pair(*it, E[*it].congestion_ratio) );
    }
  }
  //  sort ( congested_edges.begin(), congested_edges.end(), compare_pair_int_double );
  //  cout << "Congested Edges NO. = " << congested_edges.size() << endl;
}
//
//  Basic 
//
bool  CGR::isInSegment ( int vxs, Segment seg ) 
{
  if ( vxs < seg.s || vxs > seg.t )
    return  false;
  bool  isNodeInSegment;
  if ( seg.type == "HORIZ" ) {
    isNodeInSegment = (V[vxs].x>=V[seg.s].x && V[vxs].x<=V[seg.t].x && V[vxs].y==V[seg.s].y && V[vxs].z==V[seg.s].z);
  }
  if ( seg.type == "VERT" ) {
    isNodeInSegment = (V[vxs].y>=V[seg.s].y && V[vxs].y<=V[seg.t].y && V[vxs].z==V[seg.s].z && V[vxs].x==V[seg.s].x);
  }
  if ( seg.type == "VIA" ) {
    isNodeInSegment = (V[vxs].z>=V[seg.s].z && V[vxs].z<=V[seg.t].z && V[vxs].x==V[seg.s].x && V[vxs].y==V[seg.s].y);
  }
  return  isNodeInSegment;
}

int     CGR::return_edge_id ( int i, int j, int k, int type ) 
{
  if ( type==0 )
    return  k*edgeCount2D+j*(fullGridX-1)+i;
  if ( type==1 )
    return  k*edgeCount2D+ fullGridY*(fullGridX-1) + i*(fullGridY-1)+ j;
  if ( type==2 )
    return  numLayers*edgeCount2D + ( j*fullGridX + i ) * (numLayers-1) + k;
  exit(0);
}

bool    CGR::isInBox ( int node_id, GridBox &box ) 
{
  int x_coordinate = V[node_id].x;
  int y_coordinate = V[node_id].y;
  return  ( x_coordinate >= box.x0 && x_coordinate < box.x0+box.dx && y_coordinate >= box.y0 && y_coordinate < box.y0+box.dy );
}

bool    CGR::isRouted ( CoarseNet &net ) 
{
  set <int> unvisited_terminal_pins = net.terminal_pins;

  for ( vector<Segment>::iterator itr_segment = net.sg.begin(); itr_segment != net.sg.end(); itr_segment ++ ) {
    for ( set<int>::iterator itr_pin = unvisited_terminal_pins.begin(); itr_pin!= unvisited_terminal_pins.end(); itr_pin ++ ) {
      if ( isInSegment ( *itr_pin, *itr_segment ) == true )
        unvisited_terminal_pins.erase ( itr_pin );
    }
  }
  return  unvisited_terminal_pins.empty();
}

int   CGR::selfIncrementalBoundingBox ( GridBox &box ) 
{
  int lowerleft_xcoordinate = box.x0;
  int lowerleft_ycoordinate = box.y0;
  int width = box.dx;
  int height = box.dy;
  int box_util = box.util;
  int box_cap = box.cap;

  vector< pair<string,double> > added_L_congestion;
  map < string, int > L_shape_segment_added_util;
  map < string, int > L_shape_segment_added_cap;

  //  
  if ( box.x0-1 >= 0 && box.y0-1 >= 0 ) {
    int total_util = 0;
    int total_cap = 0;
    //  lower segment
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      int node_id = xyzToGridNetTerminalId ( i,box.y0-1, 0 );
      total_util += E_2d[V[node_id].ne[1]].util;
      total_util += E_2d[V[node_id].ne[2]].util;
      total_cap += E_2d[V[node_id].ne[1]].cap;
      total_cap += E_2d[V[node_id].ne[2]].cap;
    }
    //  left segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      int node_id = xyzToGridNetTerminalId ( box.x0-1, j, 0 );
      total_util += E_2d[V[node_id].ne[0]].util;
      total_util += E_2d[V[node_id].ne[3]].util;
      total_cap += E_2d[V[node_id].ne[0]].cap;
      total_cap += E_2d[V[node_id].ne[3]].cap;
    }

    added_L_congestion.push_back ( make_pair ( "lower left", (double)(box_util+total_util)/(box_cap+total_cap)) );
    L_shape_segment_added_util.insert ( make_pair ( "lower left", total_util ) );
    L_shape_segment_added_cap.insert ( make_pair ( "lower left", total_cap ) );
  }
  //
  if ( box.x0-1 >= 0 && box.y0+box.dy < fullGridY ) {
    int total_util = 0;
    int total_cap = 0;
    // top segment 
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      int node_id = xyzToGridNetTerminalId ( i, box.y0+box.dy, 0 );
      total_util += E_2d[V[node_id].ne[1]].util;
      total_util += E_2d[V[node_id].ne[3]].util;
      total_cap += E_2d[V[node_id].ne[1]].cap;
      total_cap += E_2d[V[node_id].ne[3]].cap;
    }
    //  left segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      int node_id = xyzToGridNetTerminalId ( box.x0-1, j, 0 );
      total_util += E_2d[V[node_id].ne[0]].util;
      total_util += E_2d[V[node_id].ne[2]].util;
      total_cap += E_2d[V[node_id].ne[0]].cap;
      total_cap += E_2d[V[node_id].ne[2]].cap;
    }

    added_L_congestion.push_back ( make_pair ( "top left", (double)(box_util+total_util)/(box_cap+total_cap)) );
    L_shape_segment_added_util.insert ( make_pair ( "top left", total_util ) );
    L_shape_segment_added_cap.insert ( make_pair ( "top left", total_cap ) );
  }
  //
  if ( box.x0+box.dx < fullGridX && box.y0-1 >= 0 ) {
    int total_util = 0;
    int total_cap = 0;
    //  lower segment
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      int node_id = xyzToGridNetTerminalId ( i,box.y0-1, 0 );
      assert ( V[node_id].ne[0] != -1 );
      assert ( V[node_id].ne[2] != -1 );
      total_util += E_2d[V[node_id].ne[0]].util;
      total_util += E_2d[V[node_id].ne[2]].util;
      total_cap += E_2d[V[node_id].ne[0]].cap;
      total_cap += E_2d[V[node_id].ne[2]].cap;
    }
    //  right segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      int node_id = xyzToGridNetTerminalId ( box.x0+box.dx, j, 0 );
      assert ( V[node_id].ne[1] != -1 );
      assert ( V[node_id].ne[3] != -1 );
      total_util += E_2d[V[node_id].ne[1]].util;
      total_util += E_2d[V[node_id].ne[3]].util;
      total_cap += E_2d[V[node_id].ne[1]].cap;
      total_cap += E_2d[V[node_id].ne[3]].cap;
    }

    added_L_congestion.push_back ( make_pair ( "lower right", (double)(box_util+total_util)/(box_cap+total_cap)) );
    L_shape_segment_added_util.insert ( make_pair ( "lower right", total_util ) );
    L_shape_segment_added_cap.insert ( make_pair ( "lower right", total_cap ) );
  }
  //
  if ( box.x0+box.dx < fullGridX && box.y0+box.dy < fullGridY ) {
    int total_util = 0;
    int total_cap = 0;
    //  top segment
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      int node_id = xyzToGridNetTerminalId ( i, box.y0+box.dy, 0 );
      total_util += E_2d[V[node_id].ne[0]].util;
      total_util += E_2d[V[node_id].ne[3]].util;
      total_cap += E_2d[V[node_id].ne[0]].cap;
      total_cap += E_2d[V[node_id].ne[3]].cap;
    }
    //  right segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      int node_id = xyzToGridNetTerminalId ( box.x0+box.dx, j, 0 );
      total_util += E_2d[V[node_id].ne[1]].util;
      total_util += E_2d[V[node_id].ne[2]].util;
      total_cap += E_2d[V[node_id].ne[1]].cap;
      total_cap += E_2d[V[node_id].ne[2]].cap;
    }

    added_L_congestion.push_back ( make_pair ( "top right", (double)(box_util+total_util)/(box_cap+total_cap)) );
    L_shape_segment_added_util.insert ( make_pair ( "top right", total_util ) );
    L_shape_segment_added_cap.insert ( make_pair ( "top right", total_cap ) );
  }

  //  cout << "number of qualified candidate segments = " << added_L_congestion.size() << endl;
  //  for ( vector<pair<string,double> >::iterator it = added_L_congestion.begin(); it < added_L_congestion.end(); it++ )
  //    cout << (*it).first << "\t" << (*it).second << endl;
  sort ( added_L_congestion.begin(), added_L_congestion.end(), compare_pair_string_double );
  //  cout << "Adding segment : " << added_L_congestion[0].first << endl;
  if ( added_L_congestion[0].first == "lower left" ) {
    box.util += L_shape_segment_added_util["lower left"];
    box.cap += L_shape_segment_added_cap ["lower left"];
    box.x0 --; box.y0 --;
    box.dx ++; box.dy ++;
  }
  if ( added_L_congestion[0].first == "lower right" ) {
    box.util += L_shape_segment_added_util["lower right"];
    box.cap += L_shape_segment_added_cap ["lower right"];
    box.y0 --;
    box.dx ++; box.dy ++;
  }
  if ( added_L_congestion[0].first == "top left" ) {
    box.util += L_shape_segment_added_util["top left"];
    box.cap += L_shape_segment_added_cap ["top left"];
    box.x0 --; 
    box.dx ++; box.dy ++;
  }
  if ( added_L_congestion[0].first == "top right" ) {
    box.util += L_shape_segment_added_util["top right"];
    box.cap += L_shape_segment_added_cap ["top right"];
    box.dx ++; box.dy ++;
  }
  //  cout << "num of L-shape segments = " << added_L_congestion.size() << endl;
  //  for ( vector< pair<string, double> >::iterator it = added_L_congestion.begin(); it != added_L_congestion.end(); it ++ )
  //    cout << (*it).first << "\t" << (*it).second << endl;
  return  1;
}

int   CGR::selfIncrementalBoundingBox2 ( GridBox &box ) 
{
  int lowerleft_xcoordinate = box.x0;
  int lowerleft_ycoordinate = box.y0;
  int width = box.dx;
  int height = box.dy;

  vector< pair<string,double> > added_L_congestion;
  //  
  if ( box.x0-1 >= 0 && box.y0-1 >= 0 ) {
    int sum_congestion_ratio = 0;
    //  lower segment
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      int node_id = xyzToGridNetTerminalId ( i,box.y0-1, 0 );
      sum_congestion_ratio += E_2d[V[node_id].ne[1]].congestion_ratio;
      sum_congestion_ratio += E_2d[V[node_id].ne[2]].congestion_ratio;
    }
    //  left segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      int node_id = xyzToGridNetTerminalId ( box.x0-1, j, 0 );
      sum_congestion_ratio += E_2d[V[node_id].ne[0]].congestion_ratio;
      sum_congestion_ratio += E_2d[V[node_id].ne[3]].congestion_ratio;
    }
    added_L_congestion.push_back ( make_pair ( "lower left", sum_congestion_ratio) );
  }
  //
  if ( box.x0-1 >= 0 && box.y0+box.dy < fullGridY ) {
    int sum_congestion_ratio = 0;
    // top segment 
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      int node_id = xyzToGridNetTerminalId ( i, box.y0+box.dy, 0 );
      sum_congestion_ratio += E_2d[V[node_id].ne[1]].congestion_ratio;
      sum_congestion_ratio += E_2d[V[node_id].ne[3]].congestion_ratio;
    }
    //  left segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      int node_id = xyzToGridNetTerminalId ( box.x0-1, j, 0 );
      sum_congestion_ratio += E_2d[V[node_id].ne[0]].congestion_ratio;
      sum_congestion_ratio += E_2d[V[node_id].ne[2]].congestion_ratio;
    }
    added_L_congestion.push_back ( make_pair ( "top left", sum_congestion_ratio ) );
  }
  //
  if ( box.x0+box.dx < fullGridX && box.y0-1 >= 0 ) {
    int sum_congestion_ratio = 0;
    //  lower segment
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      int node_id = xyzToGridNetTerminalId ( i,box.y0-1, 0 );
      assert ( V[node_id].ne[0] != -1 );
      assert ( V[node_id].ne[2] != -1 );
      sum_congestion_ratio += E_2d[V[node_id].ne[0]].congestion_ratio;
      sum_congestion_ratio += E_2d[V[node_id].ne[2]].congestion_ratio;
    }
    //  right segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      int node_id = xyzToGridNetTerminalId ( box.x0+box.dx, j, 0 );
      assert ( V[node_id].ne[1] != -1 );
      assert ( V[node_id].ne[3] != -1 );
      sum_congestion_ratio += E_2d[V[node_id].ne[1]].congestion_ratio;
      sum_congestion_ratio += E_2d[V[node_id].ne[3]].congestion_ratio;
    }
    added_L_congestion.push_back ( make_pair ( "lower right", sum_congestion_ratio) );
  }
  //
  if ( box.x0+box.dx < fullGridX && box.y0+box.dy < fullGridY ) {
    int sum_congestion_ratio = 0;
    //  top segment
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      int node_id = xyzToGridNetTerminalId ( i, box.y0+box.dy, 0 );
      sum_congestion_ratio += E_2d[V[node_id].ne[0]].congestion_ratio;
      sum_congestion_ratio += E_2d[V[node_id].ne[3]].congestion_ratio;
    }
    //  right segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      int node_id = xyzToGridNetTerminalId ( box.x0+box.dx, j, 0 );
      sum_congestion_ratio += E_2d[V[node_id].ne[1]].congestion_ratio;
      sum_congestion_ratio += E_2d[V[node_id].ne[2]].congestion_ratio;
    }
    added_L_congestion.push_back ( make_pair ( "top right", sum_congestion_ratio ) );
  }

  sort ( added_L_congestion.begin(), added_L_congestion.end(), compare_pair_string_double );

  if ( added_L_congestion[0].first == "lower left" ) {
    box.x0 --; box.y0 --;
    box.dx ++; box.dy ++;
  }
  if ( added_L_congestion[0].first == "lower right" ) {
    box.y0 --;
    box.dx ++; box.dy ++;
  }
  if ( added_L_congestion[0].first == "top left" ) {
    box.x0 --; 
    box.dx ++; box.dy ++;
  }
  if ( added_L_congestion[0].first == "top right" ) {
    box.dx ++; box.dy ++;
  }
  box.congestion_ratio = (double)((( box.congestion_ratio )* ( box.dx*(box.dy-1)+box.dy*(box.dx-1) ) + added_L_congestion[0].second )/( box.dx*(box.dy+1)+box.dy*(box.dx+1) ));
  return  1;
}

GridBox CGR::buildBoudingBoxFromCongestedEdge ( int edge_id ) 
{
  int srcPoint = E_2d[edge_id].srcPoint;
  int snkPoint = E_2d[edge_id].snkPoint;
  GridBox   box;
  box.set_values ( V[srcPoint].x, V[srcPoint].y, V[snkPoint].x-V[srcPoint].x+1, V[snkPoint].y-V[srcPoint].y+1 );
  box.congestion_ratio = E_2d[edge_id].congestion_ratio;
  GridBox &ref_box = box;
  for ( int itr_num = 0; itr_num < 4; itr_num++ ) {
    int status = selfIncrementalBoundingBox ( ref_box );
    if ( status == 0 ) {
      cout << "ERROR : self increment bounding box" << endl;
      exit(0);
    }
  }
  return  box;
}

GridBox CGR::buildBoudingBoxFromCongestedEdge2 ( int edge_id ) 
{ int srcPoint = E_2d[edge_id].srcPoint;
  int snkPoint = E_2d[edge_id].snkPoint;
  GridBox   box;
  int llx = max(V[srcPoint].x-1, 0);
  int lly = max(V[srcPoint].y-1, 0);
  int hrx = min(V[snkPoint].x+1, fullGridX-1);
  int hry = min(V[snkPoint].y+1, fullGridY-1);
  box.set_values ( llx, lly, hrx-llx+1, hry-lly+1 );
  return  box;
}

int     CGR::getDeadEdgeNumberOfBox (void)
{
  int number_dead_edges = 0;
  for ( vector<int>::iterator itr = E_id.begin(); itr < E_id.end(); itr++ )
    if ( E[*itr].cap==0 && E[*itr].util>0 ) {
      number_dead_edges++;
    }
  return  number_dead_edges;
}

void    CGR::countNumCongestedEdges (void)
{
  for ( vector<GridEdge>::iterator it = E.begin(); it != E.end(); it++ )
    if ( (*it).cap!=0 && (*it).util!=0 && (*it).type!="VIA" && (*it).congestion_ratio > TCR ) {
      all_the_congested_edges.insert ( it - E.begin() );
    }
}

int     CGR::getBoxes ( void ) 
{
  E_2d.clear();
  V_2d.clear();
  vector< pair<int, double> > edge_congestion_map;
  /* 2-D projection */
  for ( int j=0; j< fullGridY; j++ )
    for ( int i=0; i< fullGridX-1; i++ )    {
      GridEdge    e;
      e.id = return_edge_id ( i,j,0,0 );
      e.srcPoint = xyzToGridNetTerminalId( i,j,0 ); 
      e.snkPoint = xyzToGridNetTerminalId( i+1,j,0 ); 
      e.type = "HORIZ";
      e.congestion_ratio = 0;
      for ( int k=0; k<numLayers; k++ ) {
        int edge_id = return_edge_id ( i,j,k,0 );
        if ( E[edge_id].cap!=0 && E[edge_id].util!=0 && E[edge_id].congestion_ratio > e.congestion_ratio )
          e.congestion_ratio = E[edge_id].congestion_ratio;
      }
      edge_congestion_map.push_back ( make_pair ( e.id, e.congestion_ratio ) );
      E_2d.push_back ( e );
    }
  for ( int i=0; i<fullGridX; i++ )   
    for ( int j=0; j< fullGridY-1; j++ )  {
      GridEdge    e;
      e.id = return_edge_id ( i,j,0,1 );
      e.srcPoint = xyzToGridNetTerminalId( i,j,0 ); 
      e.snkPoint = xyzToGridNetTerminalId( i,j+1,0 ); 
      e.type = "VERT";
      e.congestion_ratio = 0;
      for ( int k=0; k<numLayers; k++ ) {
        int edge_id = return_edge_id ( i,j,k,1 );
        if ( E[edge_id].cap!=0 && E[edge_id].util!=0 && E[edge_id].congestion_ratio > e.congestion_ratio )
          e.congestion_ratio = E[edge_id].congestion_ratio;
      }
      edge_congestion_map.push_back ( make_pair ( e.id, e.congestion_ratio ) );
      E_2d.push_back ( e );
    }
  /* sort congested edges in order of their congestion ratio */
  sort ( edge_congestion_map.begin(), edge_congestion_map.end(), compare_pair_int_double );

  cout << "2-D Projection Finished" << endl;

  /* the first elements of pair in bounding_boxes_to_fix are sorted according to congestion ratio */
  bounding_boxes_to_fix.clear();

  for ( vector< pair<int,double> >::iterator it = edge_congestion_map.begin(); it < edge_congestion_map.end() && (*it).second>TCR; it ++ ) {
    int s = E_2d[(*it).first].srcPoint;
    int t = E_2d[(*it).first].snkPoint;
    GridBox box = buildBoudingBoxFromCongestedEdge ( (*it).first );
    bounding_boxes_to_fix.push_back ( make_pair((*it).first, box) );
  }

  cout << bounding_boxes_to_fix.size() << endl;
  findNonOverlappingBoxes2 ();
  return 1;
}

int     CGR::getBoxes2 ( void ) 
{
  E_2d.clear();
  V_2d.clear();
  vector< pair<int, double> > edge_congestion_map;
  /* 2-D projection */
  for ( int j=0; j< fullGridY; j++ )
    for ( int i=0; i< fullGridX-1; i++ )    {
      GridEdge    e;
      e.id = return_edge_id ( i,j,0,0 );
      e.srcPoint = xyzToGridNetTerminalId( i,j,0 ); 
      e.snkPoint = xyzToGridNetTerminalId( i+1,j,0 ); 
      e.type = "HORIZ";
      e.congestion_ratio = 0;
      for ( int k=0; k<numLayers; k++ ) {
        int edge_id = return_edge_id ( i,j,k,0 );
        if ( E[edge_id].cap>0 && e.congestion_ratio < E[edge_id].congestion_ratio )
          e.congestion_ratio = E[edge_id].congestion_ratio;
      }
      edge_congestion_map.push_back ( make_pair ( e.id, e.congestion_ratio ) );
      E_2d.push_back ( e );
    }
  for ( int i=0; i<fullGridX; i++ )   
    for ( int j=0; j< fullGridY-1; j++ )  {
      GridEdge    e;
      e.id = return_edge_id ( i,j,0,1 );
      e.srcPoint = xyzToGridNetTerminalId( i,j,0 ); 
      e.snkPoint = xyzToGridNetTerminalId( i,j+1,0 ); 
      e.type = "VERT";
      e.congestion_ratio = 0;
      for ( int k=0; k<numLayers; k++ ) {
        int edge_id = return_edge_id ( i,j,k,1 );
        if ( E[edge_id].cap>0 && e.congestion_ratio < E[edge_id].congestion_ratio )
          e.congestion_ratio = E[edge_id].congestion_ratio;
      }
      edge_congestion_map.push_back ( make_pair ( e.id, e.congestion_ratio ) );
      E_2d.push_back ( e );
    }
  /* sort congested edges in order of their congestion ratio */
  sort ( edge_congestion_map.begin(), edge_congestion_map.end(), compare_pair_int_double );

  cout << "2-D Projection Finished" << endl;

  /* the first elements of pair in bounding_boxes_to_fix are sorted according to congestion ratio */
  bounding_boxes_to_fix.clear();

  for ( vector< pair<int,double> >::iterator it = edge_congestion_map.begin(); 
      it < edge_congestion_map.begin()+1000 && (*it).second>TCR; it ++ ) {
    int s = E_2d[(*it).first].srcPoint;
    int t = E_2d[(*it).first].snkPoint;
    GridBox box = buildBoudingBoxFromCongestedEdge ( (*it).first );
    //    box.printBox();
    bounding_boxes_to_fix.push_back ( make_pair((*it).first, box) );
  }
  cout << bounding_boxes_to_fix.size() << endl;
  cout << E_2d.size() << endl;

  //  compute ACE(20%) for each box
  for ( vector<pair<int, GridBox> >::iterator it = bounding_boxes_to_fix.begin(); it != bounding_boxes_to_fix.end(); it++ ) {
    cout << it - bounding_boxes_to_fix.begin() << " " ;
    boxes_to_rank.push_back ( make_pair((*it).second, getACE20( (*it).second) ) );
    (*it).second.printBox();
    cout << getACE20 ((*it).second) << endl;
  }
  cout << "flag 2" << endl;
  //  rank all the boxes
  sort ( boxes_to_rank.begin(), boxes_to_rank.end(), compare_pair_GridBox_double );
  cout << "flag 3" << endl;
  //  print the top 20 boxes
  for ( int i=0; i < 20; i++ ) {
    boxes_to_rank[i].first.printBox();
    cout << boxes_to_rank[i].second << endl;
  }
  return 1;
}

double  CGR::getACE20 ( GridBox &box )
{
  vector<double>  congestion_ratio_of_edges;
  for ( int z=0; z < numLayers; z++ ) {
    for ( int y=box.y0; y < box.y0+box.dy; y++ )
      for ( int x=box.x0; x < box.x0+box.dx-1; x++ ){
        int  e = return_edge_id ( x, y, z, 0 );
        if ( E[e].cap>0 && E[e].util>0 )
          congestion_ratio_of_edges.push_back ( (double)(E[e].blk+E[e].util)/E[e].cap );
      }
    for ( int x=box.x0; x < box.x0+box.dx; x++ )
      for ( int y=box.y0; y < box.y0+box.dy-1; y++ ){
        int  e = return_edge_id ( x, y, z, 2 );
        if ( E[e].cap>0 && E[e].util>0 )
          congestion_ratio_of_edges.push_back ( (double)(E[e].blk+E[e].util)/E[e].cap );
      }
  }
  for ( int z=0; z < numLayers-1; z++ ) {
    for ( int y=box.y0; y < box.y0+box.dy; y++ )
      for ( int x=box.x0; x < box.x0+box.dx; x++ ){
        int  e = return_edge_id ( x, y, z, 4 );
        if ( E[e].cap>0 && E[e].util>0 )
          congestion_ratio_of_edges.push_back ( (double)(E[e].blk+E[e].util)/E[e].cap );
      }
  }
  sort ( congestion_ratio_of_edges.begin(), congestion_ratio_of_edges.end() );
  reverse ( congestion_ratio_of_edges.begin(), congestion_ratio_of_edges.end() );
  assert ( congestion_ratio_of_edges.size() > 0 );
  double  ace20 = 0;
  for ( int i=0; i < congestion_ratio_of_edges.size() * 0.2+1; i++ )
    ace20 += congestion_ratio_of_edges[i];
  return  (double)(ace20/((int)(congestion_ratio_of_edges.size()*0.2)+1));
}

vector<int> CGR::getEdgeVectorFromSegment ( Segment segment ) 
{
  vector<int> edge_vector;
  int src_point_id = segment.s;
  int snk_point_id = segment.t;
  if ( segment.type == "HORIZ" ) {
    for ( int i=V[src_point_id].x; i<V[snk_point_id].x; i++ )
      edge_vector.push_back ( return_edge_id(i, V[src_point_id].y, V[src_point_id].z, 0) );
    return  edge_vector;
  }
  if ( segment.type == "VERT" ) {
    for ( int j=V[src_point_id].y; j<V[snk_point_id].y; j++ )
      edge_vector.push_back ( return_edge_id(V[src_point_id].x, j, V[src_point_id].z, 1) );
    return  edge_vector;
  }
  if ( segment.type == "VIA" ) {
    for ( int k=V[src_point_id].z; k<V[snk_point_id].z; k++ )
      edge_vector.push_back ( return_edge_id(V[src_point_id].x, V[src_point_id].y, k, 2) );
    return  edge_vector;
  }
}

vector<int> CGR::getNodeVectorFromSegment ( Segment segment ) 
{
  vector<int> node_vector;
  int src_point_id = segment.s;
  int snk_point_id = segment.t;
  if ( segment.type == "HORIZ" ) {
    for ( int i=V[src_point_id].x; i<=V[snk_point_id].x; i++ )
      node_vector.push_back ( xyzToGridNetTerminalId ( i, V[src_point_id].y, V[src_point_id].z ) );
    return  node_vector;
  }
  if ( segment.type == "VERT" ) {
    for ( int j=V[src_point_id].y; j<=V[snk_point_id].y; j++ )
      node_vector.push_back ( xyzToGridNetTerminalId(V[src_point_id].x, j, V[src_point_id].z) );
    return  node_vector;
  }
  if ( segment.type == "VIA" ) {
    for ( int k=V[src_point_id].z; k<=V[snk_point_id].z; k++ )
      node_vector.push_back ( xyzToGridNetTerminalId(V[src_point_id].x, V[src_point_id].y, k) );
    return  node_vector;
  }
}

vector<int> CGR::getNodeVectorBetween ( int s, int t ) 
{
  vector<int> node_vector;
  int src_point_id = s;
  int snk_point_id = t;
  if ( V[s].x==V[t].x ) {
    for ( int i=V[src_point_id].x; i<=V[snk_point_id].x; i++ )
      node_vector.push_back ( xyzToGridNetTerminalId ( i, V[src_point_id].y, V[src_point_id].z ) );
    return  node_vector;
  }
  if ( V[s].y==V[t].y ) {
    for ( int j=V[src_point_id].y; j<=V[snk_point_id].y; j++ )
      node_vector.push_back ( xyzToGridNetTerminalId(V[src_point_id].x, j, V[src_point_id].z) );
    return  node_vector;
  }
  if ( V[s].z==V[t].z ) {
    for ( int k=V[src_point_id].z; k<=V[snk_point_id].z; k++ )
      node_vector.push_back ( xyzToGridNetTerminalId(V[src_point_id].x, V[src_point_id].y, k) );
    return  node_vector;
  }
}

bool    CGR::isInNet ( int edge_id_number, CoarseNet &net ) 
{
  for ( int i=0; i<net.sg.size(); i++ )
    if ( isInNetSegment ( edge_id_number, net.sg[i] ) )
      return  true;
  return  false;
}

bool    CGR::isInNetSegment ( int edge_id_number, Segment segment_in_net ) 
{
  int src_id_in_edge = E[edge_id_number].srcPoint;
  int snk_id_in_edge = E[edge_id_number].snkPoint;
  int src_id_in_segment = segment_in_net.s;
  int snk_id_in_segment = segment_in_net.t;
  if ( segment_in_net.type == "HORIZ" ) 
    return  ( V[src_id_in_edge].y==V[src_id_in_segment].y && V[src_id_in_edge].z==V[src_id_in_segment].z
        && V[src_id_in_edge].x>=V[src_id_in_segment].x && V[src_id_in_edge].x<=V[src_id_in_segment].x );
  if ( segment_in_net.type == "VERT" ) 
    return  ( V[src_id_in_edge].x==V[src_id_in_segment].x && V[src_id_in_edge].z==V[src_id_in_segment].z
        && V[src_id_in_edge].y>=V[src_id_in_segment].y && V[src_id_in_edge].y<=V[src_id_in_segment].y );
  if ( segment_in_net.type == "VIA" ) 
    return  ( V[src_id_in_edge].x==V[src_id_in_segment].x && V[src_id_in_edge].y==V[src_id_in_segment].y
        && V[src_id_in_edge].z>=V[src_id_in_segment].z && V[src_id_in_edge].z<=V[src_id_in_segment].z );
}

bool    CGR::isOverlappingBoxes ( GridBox box1, GridBox box2 ) 
{
  return ( box1.x0 <= box2.x0+box2.dx && box2.x0 <= box1.x0+box1.dx && box1.y0 <= box2.y0+box2.dy && box2.y0 <= box1.y0+box1.dy );
}

bool  CGR::isConnected ( CoarseNet &net ) 
{
  buildGraph ( net.terminal_pins, net.edges );
  return isValidRoute ( net.terminal_pins ) ;
}
//  {
//  /* build a graph */
//  map<int,set<int> >  tree_in_adjacent_form;
//  for ( set<int>::iterator itr_edge = net.edges.begin(); itr_edge != net.edges.end(); itr_edge++ ) {
//    int s = E[(*itr_edge)].srcPoint;
//    int t = E[(*itr_edge)].snkPoint;
//    if ( tree_in_adjacent_form.count ( s ) != 0 )
//      tree_in_adjacent_form[s].insert ( t );
//    else {
//      set<int>  myset;
//      myset.insert ( t );
//      tree_in_adjacent_form.insert ( make_pair(s, myset) );
//    }
//    if ( tree_in_adjacent_form.count ( t ) != 0 )
//      tree_in_adjacent_form[t].insert ( s );
//    else {
//      set<int>  myset;
//      myset.insert ( s );
//      tree_in_adjacent_form.insert ( make_pair(t, myset) );
//    }
//  }
//  /* implement a BFS and check if the number of visited nodes equals to node number of graph */
//  queue<int>  Q;
//  set<int>    V;
//  int v = (*(tree_in_adjacent_form.begin())).first;  // this can be any node in graph
//  Q.push ( v );
//  V.insert ( v );
//  while ( Q.empty()==false ) {
//    int t = Q.front();
//    Q.pop();
//    for ( set<int>::iterator it = tree_in_adjacent_form[t].begin(); it!= tree_in_adjacent_form[t].end(); it++ ) {
//      if ( V.find(*it)==V.end() ) {
//        V.insert ( *it );
//        Q.push ( *it );
//      }
//    }
//  }
//  return  V.size()==tree_in_adjacent_form.size();
//}

void  CGR::splitBoundaryTerminalPins ( CoarseNet &net ) 
{
  if ( net.boundary_terminal_pins.size() == 0 )
    return;
  /* use segment vector to determine whether they are connected or not */
  /* build a graph */
  map<int,set<int> >  tree_in_adjacent_form;
  for ( size_t i=0; i< net.sg.size(); i++ ) {
    vector<int> edge_vector = getEdgeVectorFromSegment ( net.sg[i] );
    for ( vector<int>::iterator it = edge_vector.begin(); it!= edge_vector.end(); it++ ) {
      int s = E[(*it)].srcPoint;
      int t = E[(*it)].snkPoint;
      if ( tree_in_adjacent_form.count ( s ) != 0 )
        tree_in_adjacent_form[s].insert ( t );
      else {
        set<int>  myset;
        myset.insert ( t );
        tree_in_adjacent_form.insert ( make_pair(s, myset) );
      }
      if ( tree_in_adjacent_form.count ( t ) != 0 )
        tree_in_adjacent_form[t].insert ( s );
      else {
        set<int>  myset;
        myset.insert ( s );
        tree_in_adjacent_form.insert ( make_pair(t, myset) );
      }
    }
  }
  set<int>    U;  //  unvisited boundary pins
  U.insert ( net.boundary_terminal_pins.begin(), net.boundary_terminal_pins.end() );
  while ( U.empty()==false ) {
    set<int>  connected_boundary_pins;
    int v = *(U.begin());
    set<int>    V;
    queue<int>  Q;
    Q.push ( v );
    V.insert ( v );
    while ( Q.empty()==false ) {
      int t = Q.front();
      Q.pop();
      for ( set<int>::iterator it = tree_in_adjacent_form[t].begin(); it!= tree_in_adjacent_form[t].end(); it++ ) {
        if ( V.find(*it)==V.end() ) {
          V.insert ( *it );
          Q.push ( *it );
        }
      }
    }
    for ( set<int>::iterator it = U.begin(); it != U.end(); ) {
      if ( V.find(*it) != V.end() ) {
        U.erase ( it++ );
        connected_boundary_pins.insert ( *it );
      }
      else
        ++it;
    }
    net.disjoint_boundary_terminal_pins.push_back ( connected_boundary_pins );
  }
}

bool    CGR::isLegalRoute ( CoarseNet &net ) 
{
  if ( isRouted ( net ) == false ) {
    /* all the terminal pins are visited */
    cout<<"ERROR : isRouted false "<<endl;
  }
  if ( isConnected( net ) == false ) {
    /* all the edges are connected */
    cout<<"ERROR : isConnected false "<<endl;
  }
  return  isRouted(net) && isConnected(net);
}

set<int>    CGR::getFlatSegment ( int s, int t ) 
{
  if ( s>t ) {
    int tmp = s;  s = t; t = tmp;
  }
  set<int>  edge_set;
  assert (  V[s].x==V[t].x && V[s].y==V[t].y || V[s].y==V[t].y && V[s].z==V[t].z ||  V[s].x==V[t].x && V[s].z==V[t].z ); 
  /* return integer set containing all the edges connecting two given nodes in form of flat segment */
  if ( V[s].x!=V[t].x ) {
    for ( int i = V[s].x; i < V[t].x; i++ ) 
      edge_set.insert ( return_edge_id ( i, V[s].y, V[s].z, 0 ) );
  }
  if ( V[s].y!=V[t].y ) {
    for ( int i = V[s].y; i < V[t].y; i++ ) 
      edge_set.insert ( return_edge_id ( V[s].x, i, V[s].z, 1 ) );
  }
  if ( V[s].z!=V[t].z ) {
    for ( int i = V[s].z; i < V[t].z; i++ ) 
      edge_set.insert ( return_edge_id ( V[s].x, V[s].y, i, 2 ) );
  }
  return  edge_set;
}

vector<int>   CGR::getAdjacentNodes ( vector<int> nodes, set<int> edges, int u ) 
{
  vector<int> neighbor_nodes;
  for ( int i=0; i< 6; i++ ) {
    if ( V[u].ne[i]!=-1 && edges.find(V[u].ne[i])!=edges.end() ) {
      assert ( V[u].nv[i] != -1 );
      neighbor_nodes.push_back ( V[u].nv[i] );
    }
  }
  return  neighbor_nodes;
}

int   CGR::getEdgeIdWithTwoEndingPoints ( int u, int v ) 
{
  bool  debug_mode = false;
  if ( debug_mode ) {
    printNode(u); printNode(v);
  }
  if ( u > v ) {
    int swp = u;  u = v;  v = swp;
  }
  if ( V[u].x != V[v].x ) {
    if ( V[u].ne[0]!=V[v].ne[1] )
      return  -1;
    return  V[u].ne[0];
  }
  if ( V[u].y != V[v].y ) {
    if ( V[u].ne[2]!=V[v].ne[3] )
      return  -1;
    return  V[u].ne[2];
  }
  if ( V[u].z != V[v].z ) {
    if ( V[u].ne[4]!=V[v].ne[5] )
      return  -1;
    return  V[u].ne[4];
  }
}

int   CGR::getNearestUnvisitedTerminalPin ( vector<int> nodes, map<int, set<int> > tree_in_adjacent_form, int starting_terminal_pin, set<int> visited_nodes, set<int> &path ) 
{
  //  Assigned variables :
  //    set<int> &path
  //    of course a node id as a return !
  queue<int>  q;
  set<int>    alreadly_visited;
  map<int,int>  previous_node;

  q.push ( starting_terminal_pin );
  alreadly_visited.insert ( starting_terminal_pin );
  while ( q.size() != 0 ) {
    int t = q.front();
    q.pop();
    vector<int> adjacent_nodes ( tree_in_adjacent_form[t].begin(), tree_in_adjacent_form[t].end() );
    for ( vector<int>::iterator it = adjacent_nodes.begin(); it != adjacent_nodes.end(); it ++ ) {
      if ( alreadly_visited.find(*it)==alreadly_visited.end() ) {
        q.push ( *it );
        alreadly_visited.insert ( *it );
        previous_node.insert ( make_pair (*it,t) );
        //  check if the newly visited node is target
        int current_node = *it;
        if ( find(nodes.begin(), nodes.end(), current_node)!=nodes.end()     //  is terminal pin
            &&  visited_nodes.find(current_node)==visited_nodes.end()  ) {   //  not visited yet 
          //  back-trace
          while ( current_node != starting_terminal_pin ) {
            path.insert ( getEdgeIdWithTwoEndingPoints ( current_node, previous_node[current_node] ) );
            current_node = previous_node[current_node];
          }
          return  *it;
        }
      }
    }
  }
  return  -1;
}

void    CGR::findNonOverlappingBoxes (int number_of_boxes ) 
{
  vector<pair<int,GridBox> >::iterator it = cgr.bounding_boxes_to_fix.begin();
  while ( bounding_boxes_to_solve.size() < number_of_boxes ) {
    if ( it == cgr.bounding_boxes_to_fix.end() )
      break;
    bool  valid = true;
    for ( size_t i=0; i < bounding_boxes_to_solve.size(); i++ ) {
      if ( cgr.isOverlappingBoxes ( bounding_boxes_to_solve[i], (*it).second ) ) {
        valid = false;
        it++;
        break;
      }
    }
    if ( valid==false )
      continue;
    bounding_boxes_to_solve.push_back ( (*it).second );
    (*it).second.printBox();
  }
  cout << "Select " << bounding_boxes_to_solve.size() << " bounding boxes" << endl;
}

void    CGR::findNonOverlappingBoxes2 (void) 
{
  bounding_boxes_to_solve.clear();
  for ( vector<pair<int,GridBox> >::iterator it = bounding_boxes_to_fix.begin(); it != bounding_boxes_to_fix.end(); it++ ) {
    bool  valid = true;
    for ( size_t i=0; i < bounding_boxes_to_solve.size(); i++ ) {
      if ( cgr.isOverlappingBoxes ( bounding_boxes_to_solve[i], (*it).second ) ) {
        valid = false;
        break;
      }
    }
    if ( valid==false )
      continue;
    bounding_boxes_to_solve.push_back ( (*it).second );
  }
  cout << "Select " << bounding_boxes_to_solve.size() << " non-overlapping bounding boxes" << endl;
}

//
//  SolveBox
//
double  CGR::computeSubnetEdgeWeight ( TwoTerCoarseNet &snet ) 
{
  double  subnet_weight = 0;
  for ( set<int>::iterator itr_edge = snet.edges.begin(); itr_edge != snet.edges.end(); itr_edge++ ) {
    assert ( E[*itr_edge].cap!=0 );
    if ( E[*itr_edge].type!= "VIA" ) //  Added on NOV 10
      subnet_weight += cgr.E[*itr_edge].weight;
  }
  return  subnet_weight;
}

double  CGR::computeNetRouteWeight ( NetRoute &net_route ) 
{
  double  weight = 0;
  for ( set<int>::iterator itr_edge = net_route.edges.begin(); itr_edge != net_route.edges.end(); itr_edge++ ) {
    //    assert ( E[*itr_edge].cap!=0 );
    if ( E[*itr_edge].type!="VIA" )
      weight += cgr.E[*itr_edge].weight;
  }
  return  weight;
}

bool    CGR::isLegalSubnet ( TwoTerCoarseNet &snet ) 
{
  /* 
   * GIVEN : s, t, edges
   * GOAL :  return if the given subnet is connected 
   *  */

  /* build graph */
  map<int,set<int>  >  graph;
  for ( set<int>::iterator it= snet.edges.begin(); it!= snet.edges.end(); it++ ) {
    if ( graph.count ( E[(*it)].srcPoint )!=0 )
      graph[E[(*it)].srcPoint].insert ( E[(*it)].snkPoint );
    else {
      set<int>  myset;
      myset.insert ( E[(*it)].snkPoint );
      graph.insert ( make_pair ( E[(*it)].srcPoint, myset ) );
    }
    if ( graph.count ( E[(*it)].snkPoint )!=0 )
      graph[E[(*it)].snkPoint].insert ( E[(*it)].srcPoint );
    else {
      set<int>  myset;
      myset.insert ( E[(*it)].srcPoint );
      graph.insert ( make_pair ( E[(*it)].snkPoint, myset ) );
    }
  }   //  connectivity graph is built
  /* check if exist rebundant edges */
  if ( graph[snet.s].size()!=1 || graph[snet.t].size()!=1 ) {
    cout << "SUBNET ERROR : edge rebundancy" << endl;
    return  false;
  }
  /* check if is a legal subnet */
  int curr_node = snet.s;
  set<int>  visited_nodes;
  int num_visited_edges = 0;
  while ( num_visited_edges < snet.edges.size() ) {
    /* delete adjacent nodes that are(is) already visited */
    for ( set<int>::iterator it = graph[curr_node].begin(); it!= graph[curr_node].end(); it++ ) {
      if ( visited_nodes.find(*it) != visited_nodes.end() )
        graph[curr_node].erase (*it);
    }
    if ( graph[curr_node].size()!=1 ) {
      cout << "SUBNET ERROR : illegal subnet" << endl;
      return  false;
    }
    visited_nodes.insert ( curr_node );
    curr_node = *(graph[curr_node].begin());
    ++ num_visited_edges;
  }
  assert ( curr_node == snet.t );
  return  true;
}

string  CGR::getSubnetType ( TwoTerCoarseNet &snet ) 
{
  /* determine whether the given subnet belongs to one of the pre-defined patterns */

  /* number of flat segment == number of continuous segments */
  int number_flat_segments = 1;
  vector<int> edge_vector ( snet.edges.begin(), snet.edges.end() );
  /* determine if it is a flat segment */
  for ( size_t i=0; (i+1) != edge_vector.size(); i++ ) {
    if ( edge_vector[i+1]==edge_vector[i]+1 )
      continue;
    number_flat_segments++;
  }
  if ( number_flat_segments == 1 )
    return  "flat subnet";
  if ( number_flat_segments == 2 )
    return  "L shape subnet";
  if ( number_flat_segments == 3 )
    return  "Z or U shape subnet";

  return  "other";
}

bool  CGR::isSubnetWithFlatEndingPoints ( TwoTerCoarseNet &snet ) 
{
  return  ( V[snet.s].x==V[snet.t].x && V[snet.s].y==V[snet.t].y 
      ||  V[snet.s].y==V[snet.t].y && V[snet.s].z==V[snet.t].z 
      ||  V[snet.s].z==V[snet.t].z && V[snet.s].x==V[snet.t].x  );
}

void  CGR::findCompleteUShapeCandidateSubnetSetForFlatSegment ( GridBox &box, subnetIdType subnet_id, int M2 ) 
{
  TwoTerCoarseNet &snet = bnets[ subnet_id.first ].subnets[ subnet_id.second ];
  /* all the possible candidate subnets */
  /* key = temp_candidate_subnets 
   * value = subnet weight */
  vector<pair<TwoTerCoarseNet, double> > full_pair_candidate_subnets_to_subnet_weight;
  assert ( snet.s != snet.t );
  if ( isSubnetWithFlatEndingPoints(snet) == false ) {
    //    findCompleteZShapeCandidateSubnetSetForCommonSegment ( box, subnet_id, M2 );
    return;
  }
  /* generate candidate subnets for subnet (i,j) */
  int first_edge = *(snet.edges.begin());
  /* find candidate U-shape routes in same layer 
   * since the subnet is in BFR, so only need to figure out which edges are included in new route */
  int c1,c2,c3,c4,c5,c6;
  if ( V[snet.s].y==V[snet.t].y && V[snet.s].z==V[snet.t].z  ) {
    for ( int z1 = 0; z1 < numLayers; z1++ ) {
      if ( horizCaps[z1] == 0 )
        continue;
      for ( int y = box.y0; y < box.y0+box.dy; y++ ) {
        c1 = xyzToGridNetTerminalId ( V[snet.s].x, y, z1 );
        c2 = xyzToGridNetTerminalId ( V[snet.t].x, y, z1 );
        for ( int z2 = 0; z2 < numLayers; z2++ ) {
          if ( vertCaps[z2] == 0 )
            continue;
          TwoTerCoarseNet candidate_snet;
          candidate_snet.s = snet.s;
          candidate_snet.t = snet.t;

          c3 = xyzToGridNetTerminalId ( V[snet.s].x, V[snet.s].y, z2 );
          c4 = xyzToGridNetTerminalId ( V[snet.t].x, V[snet.t].y, z2 );
          c5 = xyzToGridNetTerminalId ( V[snet.s].x, y, z2 );
          c6 = xyzToGridNetTerminalId ( V[snet.t].x, y, z2 );
          set<int>  flat_segment_1 = getFlatSegment ( snet.s, c3 );
          set<int>  flat_segment_2 = getFlatSegment ( c3, c5 );
          set<int>  flat_segment_3 = getFlatSegment ( c5, c1 );
          set<int>  flat_segment_4 = getFlatSegment ( c1, c2 );
          set<int>  flat_segment_5 = getFlatSegment ( c2, c6 );
          set<int>  flat_segment_6 = getFlatSegment ( c6, c4 );
          set<int>  flat_segment_7 = getFlatSegment ( c4, snet.t );
          candidate_snet.edges.insert ( flat_segment_1.begin(), flat_segment_1.end() );
          candidate_snet.edges.insert ( flat_segment_2.begin(), flat_segment_2.end() );
          candidate_snet.edges.insert ( flat_segment_3.begin(), flat_segment_3.end() );
          candidate_snet.edges.insert ( flat_segment_4.begin(), flat_segment_4.end() );
          candidate_snet.edges.insert ( flat_segment_5.begin(), flat_segment_5.end() );
          candidate_snet.edges.insert ( flat_segment_6.begin(), flat_segment_6.end() );
          candidate_snet.edges.insert ( flat_segment_7.begin(), flat_segment_7.end() );
          /* push pair of new candidate subnet and its weight */
          full_pair_candidate_subnets_to_subnet_weight.push_back ( make_pair ( candidate_snet, computeSubnetEdgeWeight ( candidate_snet) ) );
        }
      }
    }
  }
  if ( V[snet.s].z==V[snet.t].z && V[snet.s].x==V[snet.t].x  ) {
    for ( int z1 = 0; z1 < numLayers; z1++ ) {
      if ( vertCaps[z1] == 0 )
        continue;
      for ( int x = box.x0; x < box.x0+box.dx; x++ ) {
        c1 = xyzToGridNetTerminalId ( x, V[snet.s].y, z1 );
        c2 = xyzToGridNetTerminalId ( x, V[snet.t].y, z1 );
        for ( int z2 = 0; z2 < numLayers; z2++ ) {
          if ( horizCaps[z2] == 0 )
            continue;
          TwoTerCoarseNet candidate_snet;
          candidate_snet.s = snet.s;
          candidate_snet.t = snet.t;

          c3 = xyzToGridNetTerminalId ( V[snet.s].x, V[snet.s].y, z2 );
          c4 = xyzToGridNetTerminalId ( V[snet.t].x, V[snet.t].y, z2 );
          c5 = xyzToGridNetTerminalId ( x, V[snet.s].y, z2 );
          c6 = xyzToGridNetTerminalId ( x, V[snet.t].y, z2 );
          set<int>  flat_segment_1 = getFlatSegment ( snet.s, c3 );
          set<int>  flat_segment_2 = getFlatSegment ( c3, c5 );
          set<int>  flat_segment_3 = getFlatSegment ( c5, c1 );
          set<int>  flat_segment_4 = getFlatSegment ( c1, c2 );
          set<int>  flat_segment_5 = getFlatSegment ( c2, c6 );
          set<int>  flat_segment_6 = getFlatSegment ( c6, c4 );
          set<int>  flat_segment_7 = getFlatSegment ( c4, snet.t );
          candidate_snet.edges.insert ( flat_segment_1.begin(), flat_segment_1.end() );
          candidate_snet.edges.insert ( flat_segment_2.begin(), flat_segment_2.end() );
          candidate_snet.edges.insert ( flat_segment_3.begin(), flat_segment_3.end() );
          candidate_snet.edges.insert ( flat_segment_4.begin(), flat_segment_4.end() );
          candidate_snet.edges.insert ( flat_segment_5.begin(), flat_segment_5.end() );
          candidate_snet.edges.insert ( flat_segment_6.begin(), flat_segment_6.end() );
          candidate_snet.edges.insert ( flat_segment_7.begin(), flat_segment_7.end() );
          /* push pair of new candidate subnet and its weight */
          full_pair_candidate_subnets_to_subnet_weight.push_back ( make_pair ( candidate_snet, computeSubnetEdgeWeight ( candidate_snet) ) );
        }
      }
    }
  }

  /* Select the top M2 candidate subnets with least weight */
  if ( M2==2 ) {
    if ( full_pair_candidate_subnets_to_subnet_weight.size()==0 )
      return;
    vector<pair<TwoTerCoarseNet, double> >::iterator itr_smallest = full_pair_candidate_subnets_to_subnet_weight.begin();
    for ( vector<pair<TwoTerCoarseNet, double> >::iterator it = full_pair_candidate_subnets_to_subnet_weight.begin(); it != full_pair_candidate_subnets_to_subnet_weight.end(); it++ ) 
      if ( (*it).second < (*itr_smallest).second ) 
        itr_smallest = it ;
    bnets[subnet_id.first].candidate_subnets[subnet_id.second].push_back ( (*itr_smallest).first );

    return;
  }
  sort ( full_pair_candidate_subnets_to_subnet_weight.begin(), full_pair_candidate_subnets_to_subnet_weight.end(), compare_subnet_with_smaller_weight );

  for ( vector<pair<TwoTerCoarseNet, double> >::iterator itr = full_pair_candidate_subnets_to_subnet_weight.begin(); 
      itr < full_pair_candidate_subnets_to_subnet_weight.begin()+M2 && itr != full_pair_candidate_subnets_to_subnet_weight.end(); itr++ )
    bnets[subnet_id.first].candidate_subnets[subnet_id.second].push_back ( (*itr).first );

}

void  CGR::findCompleteZShapeCandidateSubnetSetForCommonSegment ( GridBox &box, subnetIdType subnet_id, int M2 ) 
{
  TwoTerCoarseNet &snet = bnets[ subnet_id.first ].subnets[ subnet_id.second ];
  /* all the possible candidate subnets */
  /* key = temp_candidate_subnets 
   * value = subnet weight */
  vector<pair<TwoTerCoarseNet, double> > full_pair_candidate_subnets_to_subnet_weight;
  assert ( snet.s != snet.t );

  /* generate candidate subnets for subnet (i,j) */

  int c1,c2,c3,c4,c5,c6;
  //  
  for ( int y = min(V[snet.s].y, V[snet.t].y); y <= max(V[snet.s].y, V[snet.t].y); y++ ) {
    for ( int z1=0; z1< numLayers; z1++ ) {
      if ( vertCaps[z1]==0 )
        continue;
      for ( int z2=0; z2< numLayers; z2++ ) {
        if ( horizCaps[z2]=0 )
          continue;
        for ( int z3=0; z3< numLayers; z3++ ) {
          if ( vertCaps[z3]=0 )
            continue;
          TwoTerCoarseNet candidate_snet;
          candidate_snet.s = snet.s;
          candidate_snet.t = snet.t;

          c1 = xyzToGridNetTerminalId ( V[snet.s].x, V[snet.s].y, z1 );
          c2 = xyzToGridNetTerminalId ( V[snet.s].x, y, z1 );
          c3 = xyzToGridNetTerminalId ( V[snet.s].x, y, z2 );
          c4 = xyzToGridNetTerminalId ( V[snet.t].x, y, z2 );
          c5 = xyzToGridNetTerminalId ( V[snet.t].x, y, z3 );
          c6 = xyzToGridNetTerminalId ( V[snet.t].x, V[snet.t].y, z3 );

          set<int>  flat_segment_1 = getFlatSegment ( snet.s, c1 );
          set<int>  flat_segment_2 = getFlatSegment ( c1, c2 );
          set<int>  flat_segment_3 = getFlatSegment ( c2, c3 );
          set<int>  flat_segment_4 = getFlatSegment ( c3, c4 );
          set<int>  flat_segment_5 = getFlatSegment ( c4, c5 );
          set<int>  flat_segment_6 = getFlatSegment ( c5, c6 );
          set<int>  flat_segment_7 = getFlatSegment ( c6, snet.t );
          candidate_snet.edges.insert ( flat_segment_1.begin(), flat_segment_1.end() );
          candidate_snet.edges.insert ( flat_segment_2.begin(), flat_segment_2.end() );
          candidate_snet.edges.insert ( flat_segment_3.begin(), flat_segment_3.end() );
          candidate_snet.edges.insert ( flat_segment_4.begin(), flat_segment_4.end() );
          candidate_snet.edges.insert ( flat_segment_5.begin(), flat_segment_5.end() );
          candidate_snet.edges.insert ( flat_segment_6.begin(), flat_segment_6.end() );
          candidate_snet.edges.insert ( flat_segment_7.begin(), flat_segment_7.end() );
          /* push pair of new candidate subnet and its weight */
          full_pair_candidate_subnets_to_subnet_weight.push_back ( make_pair ( candidate_snet, computeSubnetEdgeWeight ( candidate_snet) ) );
        }
      }
    }
  }
  //
  for ( int x = min(V[snet.s].x, V[snet.t].x); x <= max(V[snet.s].x, V[snet.t].x); x++ ) {
    for ( int z1=0; z1< numLayers; z1++ ) {
      if ( horizCaps[z1]==0 )
        continue;
      for ( int z2=0; z2< numLayers; z2++ ) {
        if ( vertCaps[z2]=0 )
          continue;
        for ( int z3=0; z3< numLayers; z3++ ) {
          if ( horizCaps[z3]=0 )
            continue;
          TwoTerCoarseNet candidate_snet;
          candidate_snet.s = snet.s;
          candidate_snet.t = snet.t;

          c1 = xyzToGridNetTerminalId ( V[snet.s].x, V[snet.s].y, z1 );
          c2 = xyzToGridNetTerminalId ( x, V[snet.s].y, z1 );
          c3 = xyzToGridNetTerminalId ( x, V[snet.s].y, z2 );
          c4 = xyzToGridNetTerminalId ( x, V[snet.t].y, z2 );
          c5 = xyzToGridNetTerminalId ( x, V[snet.t].y, z3 );
          c6 = xyzToGridNetTerminalId ( V[snet.t].x, V[snet.t].y, z3 );

          set<int>  flat_segment_1 = getFlatSegment ( snet.s, c1 );
          set<int>  flat_segment_2 = getFlatSegment ( c1, c2 );
          set<int>  flat_segment_3 = getFlatSegment ( c2, c3 );
          set<int>  flat_segment_4 = getFlatSegment ( c3, c4 );
          set<int>  flat_segment_5 = getFlatSegment ( c4, c5 );
          set<int>  flat_segment_6 = getFlatSegment ( c5, c6 );
          set<int>  flat_segment_7 = getFlatSegment ( c6, snet.t );
          candidate_snet.edges.insert ( flat_segment_1.begin(), flat_segment_1.end() );
          candidate_snet.edges.insert ( flat_segment_2.begin(), flat_segment_2.end() );
          candidate_snet.edges.insert ( flat_segment_3.begin(), flat_segment_3.end() );
          candidate_snet.edges.insert ( flat_segment_4.begin(), flat_segment_4.end() );
          candidate_snet.edges.insert ( flat_segment_5.begin(), flat_segment_5.end() );
          candidate_snet.edges.insert ( flat_segment_6.begin(), flat_segment_6.end() );
          candidate_snet.edges.insert ( flat_segment_7.begin(), flat_segment_7.end() );
          /* push pair of new candidate subnet and its weight */
          full_pair_candidate_subnets_to_subnet_weight.push_back ( make_pair ( candidate_snet, computeSubnetEdgeWeight ( candidate_snet) ) );
        }
      }
    }
  }
  /* Select the top M2 candidate subnets with least weight */
  if ( M2==2 ) {
    if ( full_pair_candidate_subnets_to_subnet_weight.size()==0 )
      return;
    vector<pair<TwoTerCoarseNet, double> >::iterator itr_smallest = full_pair_candidate_subnets_to_subnet_weight.begin();
    for ( vector<pair<TwoTerCoarseNet, double> >::iterator it = full_pair_candidate_subnets_to_subnet_weight.begin(); it != full_pair_candidate_subnets_to_subnet_weight.end(); it++ ) 
      if ( (*it).second < (*itr_smallest).second ) 
        itr_smallest = it ;
    bnets[subnet_id.first].candidate_subnets[subnet_id.second].push_back ( (*itr_smallest).first );

    return;
  }
  sort ( full_pair_candidate_subnets_to_subnet_weight.begin(), full_pair_candidate_subnets_to_subnet_weight.end(), compare_subnet_with_smaller_weight );

  for ( vector<pair<TwoTerCoarseNet, double> >::iterator itr = full_pair_candidate_subnets_to_subnet_weight.begin(); 
      itr < full_pair_candidate_subnets_to_subnet_weight.begin()+M2 && itr != full_pair_candidate_subnets_to_subnet_weight.end(); itr++ )
    bnets[subnet_id.first].candidate_subnets[subnet_id.second].push_back ( (*itr).first );

}

void  selfIncremental ( vector<int> &counter, vector<int> bit_range ) 
{
  assert ( counter.size()==bit_range.size() );
  if ( counter[counter.size()-1]+1 < bit_range[bit_range.size()-1] ) {
    counter[counter.size()-1]++;
    return;
  }
  /* cnt tot number */
  int total_number = 1, cnt = 0;
  for ( size_t i=0; i< bit_range.size(); i++ )
    total_number *= bit_range[i];
  /* self incremental */
  int last_bit_count = 1;
  bool  cin = true;
  while ( cnt < total_number ) {
    if ( last_bit_count > counter.size() )
      return;
    counter[counter.size()-last_bit_count] = 0;
    ++ last_bit_count;
    if ( counter[counter.size()-last_bit_count]+1 < bit_range[counter.size()-last_bit_count] ) {
      counter[counter.size()-last_bit_count]++;
      return;
    }
  }
}

vector<int>   getCandidateSubnetIdThroughCandidateNetRouteId ( vector<int> &bit_range, int net_route_id ) 
{

  vector<int> candidate_subnetid_vector;
  for ( int i = bit_range.size()-1; i >=0; i-- ) {
    candidate_subnetid_vector.push_back ( net_route_id % bit_range[i] );
    net_route_id = net_route_id / bit_range[i];
  }
  reverse ( candidate_subnetid_vector.begin(), candidate_subnetid_vector.end() );
  return  candidate_subnetid_vector;

}

void  CGR::addNewCandidateRoutes ( int bnet_index, int upper_bound_number ) 
{
  bool  debug_mode = false;
  int total_number = 1;
  /* simulate a counter */
  vector< int > selected_subnets ( bnets[bnet_index].candidate_subnets.size(), 0 );
  vector< int > &ref_selected_subnets = selected_subnets;
  vector< int > subnets_size ( bnets[bnet_index].candidate_subnets.size(), 0 );
  for ( size_t i=0; i < bnets[bnet_index].candidate_subnets.size(); i++ ) {
    subnets_size[i] = bnets[bnet_index].candidate_subnets[i].size();
    total_number *= subnets_size[i];
  }
  assert ( total_number != 0 );
  //  cout << "Before : " << bnets[bnet_index].candidate_routes.size() << endl;
  vector<pair<NetRoute, double> > route_weight;
  for ( int i=0; i < total_number; i++ ) {
    /* build candidate route with candidate subnets in index vector */
    vector<TwoTerCoarseNet> subnets;
    for ( size_t j= 0; j < bnets[bnet_index].candidate_subnets.size(); j++ )
      subnets.push_back ( bnets[bnet_index].candidate_subnets[j][ref_selected_subnets[j]] );
    /* construct a new candidate net route using subnet vector */
    NetRoute  new_netroute ( subnets );
    route_weight.push_back ( make_pair ( new_netroute, computeNetRouteWeight (new_netroute) ) );
    /* For Debug */
    if ( debug_mode ) {
      cout << "candidate net route " << i << " : edge number = "  << new_netroute.edges.size() << endl;
      cout << "number of subnets = " << bnets[bnet_index].candidate_subnets.size() << endl;
      for ( size_t j= 0; j < bnets[bnet_index].candidate_subnets.size(); j++ )
        cout << ref_selected_subnets[j] << ",";
      cout << endl;
    }
    /* End */
    selfIncremental ( ref_selected_subnets, subnets_size );
  }
  sort ( route_weight.begin(), route_weight.end(), compare_netroute_with_smaller_weight );
  for ( size_t ii=0; ii < min(total_number, upper_bound_number); ii++ )
    bnets[bnet_index].candidate_routes.push_back ( route_weight[ii].first );
  //  cout << "After : " << bnets[bnet_index].candidate_routes.size() << endl;
}

/* 
 * This function has only one parameter of a box reference and only needs to run once per box iteration.
 * Here is the function of this module : 
 *  1. Clear everything and build a vector to store all the edges in box.
 *  2. Check for all the nets and select the related nets into a vector called bnets, which stores all the modified nets. And these nets are the objectives to deal with in this box iteration. Each bnet has its own old index, terminal pin vector, segments vector and edge set, like what a common net does.
 *  */
void    CGR::genModifiedNetsInBox (GridBox &box) 
{
  /* put all the edges inside bounding box into box.E_id */
  E_id.clear();
  bool  debug_mode = false;
  for ( int k = 0; k < numLayers; k++ )   {
    for ( int j=box.y0; j< box.y0+box.dy; j++ )
      for ( int i=box.x0; i< box.x0+box.dx-1; i++ ) {
        int edge_id = return_edge_id ( i,j,k,0 );
        map_from_edge_id_to_index_in_E_id.insert (  make_pair ( edge_id, E_id.size() ) );
        E_id.push_back ( edge_id );
      }
    for ( int i=box.x0; i< box.x0+box.dx; i++ ) 
      for ( int j=box.y0; j< box.y0+box.dy-1; j++ )  {
        int edge_id = return_edge_id ( i,j,k,1 );
        map_from_edge_id_to_index_in_E_id.insert (  make_pair ( edge_id, E_id.size() ) );
        E_id.push_back ( edge_id );
      }
  }
  for ( int j=box.y0; j< box.y0+box.dy; j++ )
    for ( int i=box.x0; i< box.x0+box.dx; i++ ) 
      for ( int k = 0; k < numLayers-1; k++ )  {
        int edge_id = return_edge_id ( i,j,k,2 );
        map_from_edge_id_to_index_in_E_id.insert (  make_pair ( edge_id, E_id.size() ) );
        E_id.push_back ( edge_id );
      }
  cout << "Built Grid Box " << endl;
  /* end */
  bnets.clear();
  for ( size_t i=0; i<nets.size(); i++ ) {
    CoarseNet &ref_net = nets[i];
    CoarseNet   modified_net_inbox = genSingleModifiedNetInBox (ref_net, box );
    if ( modified_net_inbox.edges.size() > 1 ) 
      bnets.push_back ( modified_net_inbox );
  }
  cout << "In this box iteration, the number of modified nets = " << bnets.size() << endl;
  //  merge short segments into longer ones 
  for ( size_t i=0; i<bnets.size(); i++ ) {
    CoarseNet &net = bnets[i];
    vector<int>  edges ( net.edges.begin(), net.edges.end() );
    net.sg.clear();
    // the flat segment has a very good property that all the edge id numbers are continuous. 
    vector<int>::iterator itr_edge = edges.begin();
    vector<int>::iterator itr_segment_begin = edges.begin(); 
    vector<int>::iterator itr_segment_end;

    while ( (itr_edge+1) != edges.end() ) {
      if ( *(itr_edge+1) == (*itr_edge)+1 && E[*itr_edge].snkPoint==E[*(itr_edge+1)].srcPoint && E[*itr_edge].type==E[*(itr_edge+1)].type ) {
        itr_edge++;
        continue;
      }
      itr_segment_end = itr_edge;

      Segment merged_segment;
      merged_segment.s = E[*itr_segment_begin].srcPoint;
      merged_segment.t = E[*itr_segment_end].snkPoint;
      merged_segment.type = E[*itr_segment_begin].type;

      net.sg.push_back ( merged_segment );
      itr_edge++;
      itr_segment_begin = itr_edge;
    }
    itr_segment_end = itr_edge;

    Segment merged_segment;
    merged_segment.s = E[*itr_segment_begin].srcPoint;
    merged_segment.t = E[*itr_segment_end].snkPoint;
    merged_segment.type = E[*itr_segment_begin].type;

    net.sg.push_back ( merged_segment );
  }
  cout << "Segment Merged" << endl;

  for ( size_t i=0; i < bnets.size(); i++ )
    if ( debug_mode== true ) {
      cout << bnets.size() << " " << bnets[i].old_index << " " << bnets[i].terminal_pins.size() << " " << bnets[i].edges.size() << " " << isConnected(bnets[i] ) << endl;
      //      printNet ( bnets[i] );
    }
}

/* 
 * This function exectues right after "genModifiedNetsInBox" and only runs once per box iteration. This function implements the initialization for followed ILP iterations and enables the followed functions can be executed in a normal way.
 *  1. "modified_net_type" is defined for each modified net to show whether it is connected or not. Since this property will keep unchanged during the same box iteration, this vector only needs to compute once and be mapped to each modified net until the box iterations ends. And the modified net type is determined by calling function "isConnected" by reading the edge set of bnet.
 *  2. Implement net decomposition for connected modified net. Net decomposition generates a vector of subnets, each of which contains two terminal pins and edge set. 
 *  3. As for disconnected modified net, split the boundary terminal pins into disjoint sets and hold it. 
 *  4. Store base route for both types of modified nets. For type 0, the candidate route vector is initialized with a net route containing both subnet vector and edge set. Candidate subnet is also initialized with correct number of entries, each of which is a vector of one subnet, i.e. the original subnet. For type 1, there is no candidate subnets and the only item in candidate route has only a edge set.
 * */
void  CGR::initializeLoop ( ) 
{
  /* build an integer array store the net type of modified nets */

  for ( size_t i=0; i < bnets.size(); i++ ) {
    CoarseNet &net = bnets[i];

    if ( isConnected ( net ) ) {
      modified_net_type.push_back ( 0 );
      netDecomposition ( net );
      /* initialize with base route */
      bnets[i].candidate_subnets.clear();
      bnets[i].candidate_routes.clear();
      for ( size_t j = 0; j < bnets[i].subnets.size(); j++ ) {
        vector < TwoTerCoarseNet > myset;
        myset.push_back ( bnets[i].subnets[j] );
        bnets[i].candidate_subnets.push_back ( myset );
      }
      NetRoute  base_route ( bnets[i].subnets );
      //  1st element for candidate routes
      bnets[i].candidate_routes.push_back ( base_route );
    }
    else {
      modified_net_type.push_back ( 1 );
      splitBoundaryTerminalPins ( net ); 
      /* initialize with base route */
      bnets[i].candidate_routes.clear();
      NetRoute base_route;
      base_route.edges = bnets[i].edges;
      //  1st element for candidate routes
      bnets[i].candidate_routes.push_back ( base_route );
    }
  }
}

/* 
 * This function is executed right after each time ILP is solved. In previous function of solving ILP, the primal solution and dual solution will give which among candidate routes per net is selected and which value should be assigned to edge weight. 
 *  1.  re-initialize the candidate_routes with new solution
 *  2.  re-calculate congestion edges
 *  3.  re-decompose net and initialize the candidate_subnets vector
 * */
void  CGR::readyForNextIteration (GridBox &box ) 
{
  /* re-initialize the candidate routes with new solution 
   * replace the set<int> edges in net with the edge set of selected candidated route */
  for ( size_t i=0; i< bnets.size(); i++ ) {
    /* update the subnets of given net */
    bnets[i].subnets.clear();
    bnets[i].candidate_subnets.clear();
    if ( modified_net_type[i] == 0 ) {
      netDecomposition ( bnets[i] );
      /* initialize the candidate subnet */
      for ( size_t j=0; j < bnets[i].subnets.size(); j++ ) {
        TwoTerCoarseNet base_subnet;
        base_subnet.s = bnets[i].subnets[j].s;
        base_subnet.t = bnets[i].subnets[j].t;
        base_subnet.edges = bnets[i].subnets[j].edges;

        vector< TwoTerCoarseNet > candidate_snet;
        candidate_snet.push_back ( base_subnet );
        bnets[i].candidate_subnets.push_back ( candidate_snet );
      }
    }
  }

  /* re-calculate congested edges inside box */
  for ( vector<int>::iterator itr_edge = E_id.begin(); itr_edge != E_id.end(); itr_edge++ ) 
    E[*itr_edge].util = 0;
  for ( size_t i=0; i< bnets.size(); i++ ) {
    for ( set<int>::iterator itr_edge = bnets[i].edges.begin(); itr_edge!= bnets[i].edges.end(); itr_edge++ ) {
      ++ E[*itr_edge].util;
    }
  }
  for ( vector<int>::iterator itr_edge = E_id.begin(); itr_edge != E_id.end(); itr_edge++ ) 
    E[*itr_edge].congestion_ratio = (double)((E[*itr_edge].blk+E[*itr_edge].util)/(E[*itr_edge].cap+0.0001));

  findCongestedEdges ();

  subnets_to_expand.clear();
  //  for ( vector<pair<int,double> >::iterator itr_edge = congested_edges.begin(); itr_edge!=congested_edges.end(); itr_edge++ ) {
  //  APPROXIMATION 1
  for ( vector<pair<int,double> >::iterator itr_edge = congested_edges.begin(); itr_edge!=congested_edges.begin()+congested_edges.size()*A1; itr_edge++ ) {
    int edge_id = (*itr_edge).first;
    //  Added on NOV 12
    //  Only generate candidate subnets for portion of all the subnets passing through HUT edge
    int number_considered_subnets = 0;
    //  APPROXIMATION 2
    while ( number_considered_subnets < E[edge_id].util * A2 ) {
      for ( size_t i=0; i< bnets.size(); i++ ) {
        if ( modified_net_type[i] == 1 )
          continue;
        for ( size_t j=0; j< bnets[i].subnets.size(); j++ ) {
          /* 
           * for each subnet (i,j), check if it contains the edge
           */
          TwoTerCoarseNet &snet = bnets[i].subnets[j];
          bool is_edge_contained_in_this_subnet = (snet.edges.find( (*itr_edge).first ) != snet.edges.end() );
          if ( is_edge_contained_in_this_subnet == false )  //  if the subnet does not contain that edge, then nothing happens.
            continue;
          subnets_to_expand.insert ( make_pair(i,j) );
        }
      }

      ++ number_considered_subnets;
    }
  }
  int M2 = 2;
  int M3 = 4;
  for ( set<subnetIdType>::iterator itr = subnets_to_expand.begin(); itr != subnets_to_expand.end(); itr++ ) {
    /* find up to M2 candidate routes for given subnet */
    findCompleteUShapeCandidateSubnetSetForFlatSegment ( box, *itr, M2 ); 
  }
  cout << "generated at most " << M2 << " candidate routes for each subnet " << endl;
  /* 
   * This is a module connected to another module to add columns.
   */
  for ( vector<CoarseNet>::iterator itr_net = bnets.begin(); itr_net != bnets.end(); itr_net ++ ) {
    int num_routes = 1;
    int ii = itr_net - bnets.begin();
    if ( modified_net_type[ii] != 0 ) 
      continue; 
    addNewCandidateRoutes ( ii, M3 );
    /* 
       cout << "Bnet " << itr_net-bnets.begin() << endl;
       for ( vector<NetRoute>::iterator itr = (*itr_net).candidate_routes.begin(); itr< (*itr_net).candidate_routes.end(); itr++ )
       cout << "\tcandidate rout " << itr-(*itr_net).candidate_routes.begin()  << "\t" << computeNetRouteWeight ( *itr) << endl;
       */
  }

  cout << "combined at most " << M3 << " candidate net routes" << endl;
}

/* 
 *  This function merges the net edges inside bounding box together with the parts that lay outside the bounding box. Also refine the segments and check the legality. 
 * */
/* 
 *  This function only needs the edge set of the given net and generate a segment vector, each of which are made as long as possible.
 *  */
void  CGR::generateLongestSegments ( int index_in_nets ) 
{
  CoarseNet &net = nets[index_in_nets];
  vector<int>  edges ( net.edges.begin(), net.edges.end() );
  /* the flat segment has a very good property that all the edge id numbers are continuous. */
  vector<int>::iterator itr = edges.begin();
  vector<int>::iterator itr_begin = edges.begin(); 

  net.sg.clear();
  while ( itr_begin != edges.end() ) {
    if ( (itr+1) == edges.end() || (*itr)+1!=(*(itr+1)) || E[*itr].snkPoint != E[*(itr+1)].srcPoint || E[*itr].type!=E[*(itr+1)].type ) {
      assert ( E[*itr].type == E[*itr_begin].type );
      Segment merged_segment;
      merged_segment.s = E[*itr_begin].srcPoint;
      merged_segment.t = E[*itr].snkPoint;
      merged_segment.type = E[*itr].type;
      net.sg.push_back ( merged_segment );
      itr_begin = itr+1;
    }
    ++ itr;
  }
  //  printSegmentEdges ( new_segment );
}
//  
//  Performance Metric
//
/* 
 * ACE ( y% ) 
 *  the average congestion ratio of the top y% congested edges
 *  FYI : here we skip all the complete blockages, i.e. edges with zero capacity 
 *  */
double  CGR::getACE ( double y ) 
{
  vector<double>  v;
  for ( vector<GridEdge>::iterator itr_edge = E.begin(); itr_edge != E.end(); itr_edge++ ) 
    if ( (*itr_edge).cap != 0 && (*itr_edge).type!="VIA" )
      v.push_back ( (*itr_edge).congestion_ratio );
  sort ( v.begin(), v.end() );
  reverse ( v.begin(), v.end() );

  double  sum_congestion_edge = 0;
  int edge_count = 0;
  for ( size_t i=0; i < v.size() * y ; i++ ) {
    sum_congestion_edge += v[i];
    ++ edge_count;
  }
  return  (double)sum_congestion_edge / edge_count;
}

double  CGR::getNetCongestioRatio ( CoarseNet &net ) 
{
  double  g_net= 0;
  for ( set<int>::iterator itr_edge = net.edges.begin(); itr_edge != net.edges.end();  itr_edge ++ ) {
    assert ( E[*itr_edge].cap != 0 );
    if ( E[*itr_edge].type!="VIA" && E[*itr_edge].congestion_ratio > g_net )
      g_net = E[*itr_edge].congestion_ratio;
  }
  return  g_net;
}

/* 
 * WCI ( TCR ) 
 *  the number of nets with congestion ratio greater than TCR 
 *   where the congestion ratio of net is defined as the highest congestion ratio of net edges
 *    FYI : ignore congestion ratio of all the edges with zero capacity */
int CGR::getWCI ( double TCR ) 
{
  int counter = 0;
  for ( size_t i = 0; i < nets.size(); i++ ) {
    double  net_congestion_ratio = getNetCongestioRatio ( nets[i] );
    if ( net_congestion_ratio > TCR )
      counter++;
  }
  return  counter;
}

double  CGR::getGlobalObjVal (void) 
{
  double  global_obj_val = 0;
  for ( size_t i=0; i < E.size(); i++ ) {
    if ( E[i].cap==0 || E[i].type=="VIA" )
      continue;
    global_obj_val += getEdgePenalty ( i );
//    if ((double)(E[i].blk+E[i].util)/(E[i].cap+0.0001)>TCR)
//      global_obj_val += slope* ((double)(E[i].blk+E[i].util)/(E[i].cap+0.0001)-TCR);
  }
  return  global_obj_val;
}

void  CGR::getInputGlobalRoutingInfo ( std::ofstream &stream_file ) 
{
  total_overflow = 0;
  for ( size_t i = 0; i< E.size(); i++ ) {
    assert ( E[i].util == E[i].nets_passing_through.size() );
    total_overflow += (E[i].util+E[i].blk>E[i].cap)?(E[i].util+E[i].blk-E[i].cap):0 ;
  }

  stream_file << total_overflow << " ";
  stream_file << getGlobalObjVal() << " ";
  stream_file << getACE ( 0.005 ) << " ";
  stream_file << getACE ( 0.01 ) << " ";
  stream_file << getACE ( 0.02 ) << " ";
  stream_file << getACE ( 0.05 ) << " ";
  stream_file << getACE ( 0.1 ) << " ";
  stream_file << getACE ( 0.2 ) << " ";
  stream_file << getTWL () << " ";

  stream_file << endl;
}

int   CGR::countBadEdges ( void ) 
{
  /* return the number of bad edges in global grid */
  int number_bad_edges = 0;
  for ( vector<GridEdge>::iterator itr_edge = E.begin(); itr_edge != E.end(); itr_edge ++ ) {
    if ( (*itr_edge).cap==0 && (*itr_edge).util>0 )
      number_bad_edges++;
  }
  return  number_bad_edges;
}

int   CGR::countBlockedEdges (void) 
{
  /* return the number of blocked edges in global grid */
  int number_blocked_edges = 0;
  for ( vector<GridEdge>::iterator itr_edge = E.begin(); itr_edge != E.end(); itr_edge ++ ) {
    if ( (*itr_edge).blk > 0 )
      number_blocked_edges++;
  }
  return  number_blocked_edges;
}

bool  CGR::isLegalSegment ( Segment seg ) 
{
  if ( seg.type == "VIA" )
    return  true;
  if ( seg.type == "HORIZ" ) {
    return  ( horizCaps[V[seg.s].z] != 0 );
  }
  if ( seg.type == "VERT" ) {
    return  ( vertCaps[V[seg.s].z] != 0 );
  }
  cout << "ERROR" << endl;
  exit(0);
}

/*
 * this function removes given edge from given net
 */
void  CGR::releaseDeadEdgeFromNet ( int edge_id, int net_id ) 
{
  assert ( E[edge_id].type != "VIA" );
  assert ( nets[net_id].edges.find(edge_id) != nets[net_id].edges.end() ); 
  nets[net_id].edges.erase ( edge_id );
  //
  //  cout << "\tfrom " ; 
  //  printEdge ( edge_id );
  //
  E[edge_id].util --;
  E[edge_id].congestion_ratio = E[edge_id].util / ( E[edge_id].cap-E[edge_id].blk+0.0001 );
  E[edge_id].nets_passing_through.erase ( net_id );
  //  select edge with smallest congestion ratio
  int new_e = edge_id;
  for ( int e = edge_id%edgeCount2D; e < edgeCount2D * numLayers; e += edgeCount2D ) 
    if ( E[e].congestion_ratio <= E[new_e].congestion_ratio ) 
      new_e = e;
  cout << "FROM : ";  printEdge ( edge_id );
  cout << "TO : " ;   printEdge ( new_e );

  return;

  assert ( V[E[new_e].srcPoint].x==V[E[edge_id].srcPoint].x );
  assert ( V[E[new_e].snkPoint].x==V[E[edge_id].snkPoint].x );
  assert ( V[E[new_e].srcPoint].y==V[E[edge_id].srcPoint].y );
  assert ( V[E[new_e].snkPoint].y==V[E[edge_id].snkPoint].y );
  //  re-connect
  set<int>  flat_sg1 = getFlatSegment ( E[edge_id].srcPoint, E[edge_id].srcPoint );
  set<int>  flat_sg2 = getFlatSegment ( E[edge_id].snkPoint, E[edge_id].snkPoint );
  nets[net_id].edges.insert ( flat_sg1.begin(), flat_sg1.end() );
  nets[net_id].edges.insert ( new_e );
  for ( set<int>::iterator it1 = flat_sg1.begin(); it1!= flat_sg1.end(); it1++ )
    if ( nets[net_id].edges.find(*it1) == nets[net_id].edges.end() ) {
      E[*it1].util ++;
      E[*it1].congestion_ratio = (E[*it1].blk+E[*it1].util)/(E[*it1].cap+0.0001 );
      E[*it1].nets_passing_through.insert ( net_id );
    }
  for ( set<int>::iterator it2 = flat_sg1.begin(); it2!= flat_sg1.end(); it2++ )
    if ( nets[net_id].edges.find(*it2) == nets[net_id].edges.end() ) {
      E[*it2].util ++;
      E[*it2].congestion_ratio = (E[*it2].blk+E[*it2].util)/(E[*it2].cap+0.0001 );
      E[*it2].nets_passing_through.insert ( net_id );
    }
  nets[net_id].edges.insert ( flat_sg2.begin(), flat_sg2.end() );
  if ( nets[net_id].edges.find(new_e) == nets[net_id].edges.end() ) {
    E[new_e].util ++;
    E[new_e].congestion_ratio = E[new_e].util/(E[new_e].cap-E[new_e].blk+0.0001 );
    E[new_e].nets_passing_through.insert ( net_id );
  }
  //
  //  cout << "\tto " ;
  //  printEdge ( new_e );
  //
}

int   CGR::clearDeadEdges (void)
{
  int cnt = 0;
  for ( vector<GridEdge>::iterator it = E.begin(); it < E.end(); it++ )
    if ( (*it).cap==0 && (*it).util>0 ) {
      //      cout << cnt << " Edge " << (*it).id << endl;
      //      printEdge ( (*it).id );
      for ( set<int>::iterator it2 = (*it).nets_passing_through.begin(); it2 != (*it).nets_passing_through.end(); it2++ ) {
        //        cout << "\tnets passing through : " << *it2 << endl;
        releaseDeadEdgeFromNet ( (*it).id, *it2 );
      }
      cnt++;
    }
  cout << "cnt = " << cnt << endl;
  exit(0);
  return  1;
}
//
//  Echo
//
void  CGR::printSegmentEdges ( Segment segment ) 
{
  cout << "\t(" << V[segment.s].x << "," << V[segment.s].y << "," << V[segment.s].z << ") - (";
  cout << V[segment.t].x << "," << V[segment.t].y << "," << V[segment.t].z << ")" << endl;
}

void  CGR::printEdgeSet ( set<int> edges ) 
{
  for ( set<int>::iterator it = edges.begin(); it != edges.end(); it++ ) {
    int s = E[*it].srcPoint;
    int t = E[*it].snkPoint;
    printEdge ( *it );
  }
  cout << endl;
}

void  CGR::printSegmentVector ( vector<Segment> sg ) 
{
  for ( size_t  i=0;  i<sg.size();  i++ )
    printSegmentEdges ( sg[i] );
}

void  CGR::printCrossingPointsOnSegments ( CoarseNet &net ) 
{
  for ( vector<set<int> >::iterator it = net.crossing_points_on_segments.begin(); it != net.crossing_points_on_segments.end(); it++ ) {
    cout << "---" << endl;
    printSegmentEdges ( net.sg[it-net.crossing_points_on_segments.begin()] );
    for ( set<int>::iterator it2 = (*it).begin(); it2 != (*it).end(); it2++ ) {
      cout << "\t(" << V[*it2].x << "," << V[*it2].y << "," << V[*it2].z << ")" << endl;
    }
  }
}

void  CGR::printVertices ( vector<int>  vxs ) 
{
  for ( size_t i=0; i< vxs.size(); i++ )
    cout << "\t(" << V[vxs[i]].x << "," << V[vxs[i]].y << "," << V[vxs[i]].z << ")" << endl;
}

void  CGR::printEdge ( int edge_id ) 
{
  int s = E[edge_id].srcPoint;
  int t = E[edge_id].snkPoint;
  cout << "edge id : " << edge_id << "\t(" ;
  cout << V[s].x <<","<<V[s].y<<","<<V[s].z<<")-(";
  cout << V[t].x <<","<<V[t].y<<","<<V[t].z<<") : " 
    << E[edge_id].type << " " << E[edge_id].cap <<" " 
       << E[edge_id].blk << " " << E[edge_id].util << " " 
       << "\t" << E[edge_id].congestion_ratio << " " << E[edge_id].weight << endl;
}

void  CGR::printEdge_2d ( int edge_id ) 
{
  int s = E_2d[edge_id].srcPoint;
  int t = E_2d[edge_id].snkPoint;
  cout << "(" << V[s].x <<","<<V[s].y<<","<<V[s].z<<")-(";
  cout << V[t].x <<","<<V[t].y<<","<<V[t].z<<")\tcap=" << E_2d[edge_id].cap <<" util=" << E_2d[edge_id].util << "\t" << E_2d[edge_id].congestion_ratio << endl;
}

void  CGR::printNode ( int node_id ) 
{
  cout << "node : (" << V[node_id].x<<","<< V[node_id].y<<"," << V[node_id].z << ")" << endl;
}

void  CGR::printNode_2d ( int node_id ) 
{
  cout << "node : (" << V[node_id].x<<","<< V[node_id].y<< ")\t" << V_2d[node_id].util << " " << V_2d[node_id].cap << " " << V_2d[node_id].congestion_ratio << endl;
  //  for ( size_t j=0; j<4; j++ )
  //    if ( V[node_id].ne[j] != -1 )
  //      printEdge_2d ( V_2d[node_id].ne[j] );
}

void  CGR::printEdgeWeightInSegment ( Segment segment ) 
{
  vector<int> edge_vector = getEdgeVectorFromSegment ( segment );
  for ( size_t i=0; i< edge_vector.size(); i++ ) {
    int s = E[edge_vector[i]].srcPoint;
    int t = E[edge_vector[i]].snkPoint;
    cout << "(" << V[s].x <<","<<V[s].y<<","<<V[s].z<<")-(";
    cout << V[t].x <<","<<V[t].y<<","<<V[t].z<<")\t" 
      << E[edge_vector[i]].cap  << " " 
      << E[edge_vector[i]].util << " " 
      << E[edge_vector[i]].weight << endl;
  }
}

void  CGR::printIntVector ( vector<int> int_vector ) 
{
  for ( size_t i=0; i<int_vector.size(); i++ )
    cout << int_vector[i] << " ";
  cout << endl;
}

void  CGR::printNet ( CoarseNet &net ) 
{
  cout << "net id " << net.old_index << endl;
  cout << "\tterminal pins : " << endl;
  vector<int>   terminal_pins ( net.terminal_pins.begin(), net.terminal_pins.end() );
  printVertices ( terminal_pins );
  cout << "\tsegments : " << endl;
  printSegmentVector ( net.sg );
}

void  CGR::printSubNet ( TwoTerCoarseNet &snet ) 
{
  cout << "FROM : " ; printNode ( snet.s );
  cout << "TO : " ; printNode ( snet.t );
  printEdgeSet ( snet.edges );
}
//  Added on NOV6
int   CGR::getNetHPWL ( CoarseNet &net ) 
{

  int xmin = fullGridX;
  int ymin = fullGridY;
  int xmax = 0;
  int ymax = 0;

  assert ( net.terminal_pins.size()>1 );
  for ( set<int>::iterator it = net.terminal_pins.begin(); it != net.terminal_pins.end(); it++ ) {
    int x = V[*it].x;
    int y = V[*it].y;
    xmin = min(x, xmin);  xmax = max(x, xmax);
    ymin = min(y, ymin);  ymax = max(y, ymax);
  }
  assert ( xmin<=xmax && ymin<=ymax );
  return  xmax-xmin+ymax-ymin;
}

void  CGR::printNetHPWLDistribution (void)
{
  vector<int> net_hpwl;
  for ( vector<CoarseNet>::iterator it = nets.begin(); it != nets.end(); it++ )
    if ( getNetHPWL(*it)>10 )
      (*it).edges.clear();
}

set<int>  CGR::getBestRouteForFlatNet (int s, int t)
{
  set<int>  edges;
  assert ( V[s].x!=V[t].x || V[s].y!=V[t].y );
  //
  //  if horizontal
  if ( V[s].y==V[t].y ) {
    int curr_layer = V[s].z;
    //  for each non-via edge
    for ( int x = V[s].x; x < V[t].x; x++ ) {
      int z = 0;
      double  min_weight = 1000000;
      //  determine the layer of this non-via edge
      for ( int k=0; k < numLayers; k++ ) {
        if ( horizCaps[k] == 0 )
          continue;
        if ( E[return_edge_id(x, V[s].y, k, 0)].weight < min_weight ) {
          min_weight = E[return_edge_id(x,V[s].y,k,0)].weight;
          z = k;
        }
      }
      assert ( z!= 0 );
      //  insert edge
      edges.insert ( return_edge_id(x,V[s].y,z,0) );
      //  insert vias
      for ( int kk = min(curr_layer,z ); kk < max(curr_layer,z); kk++ )
        edges.insert ( return_edge_id( x, V[s].y, kk, 2 ) ); 
      curr_layer = z;
    }
    //  insert vias
    for ( int kk = min(curr_layer,V[t].z ); kk < max(curr_layer,V[t].z); kk++ )
      edges.insert ( return_edge_id( V[t].x, V[t].y, kk, 2 ) ); 
    return  edges;
  }
  //  
  //  if vertical
  if ( V[s].x==V[t].x ) {
    int curr_layer = V[s].z;
    //  for each non-via edge
    for ( int y = V[s].y; y < V[t].y; y++ ) {
      int z = 0;
      double  min_weight = 1000000;
      //  determine the layer of this non-via edge
      for ( int k=0; k < numLayers; k++ ) {
        if ( vertCaps[k] == 0 )
          continue;
        if ( E[return_edge_id(V[s].x, y, k, 1)].weight < min_weight ) {
          min_weight = E[return_edge_id(V[s].x, y, k, 1)].weight;
          z = k;
        }
      }
      assert ( z!= 0 );
      //  insert edge
      edges.insert ( return_edge_id( V[s].x, y, z, 1) );
      //  insert vias
      for ( int kk = min(curr_layer,z ); kk < max(curr_layer,z); kk++ )
        edges.insert ( return_edge_id( V[s].x, y, kk, 2 ) ); 
      curr_layer = z;
    }
    //  insert vias
    for ( int kk = min(curr_layer,V[t].z ); kk < max(curr_layer,V[t].z); kk++ )
      edges.insert ( return_edge_id( V[t].x, V[t].y, kk, 2 ) ); 
    return  edges;
  }
  cout << "ERROR" << endl;
  exit(0);
}
//
//  Added on NOV7
void    CGR::genModifiedNetsInBox2 (GridBox &box) 
{
  /* put all the edges inside bounding box into box.E_id */
  E_id.clear();
  bool  debug_mode = false;
  for ( int k = 0; k < numLayers; k++ )   {
    for ( int j=box.y0; j< box.y0+box.dy; j++ )
      for ( int i=box.x0; i< box.x0+box.dx-1; i++ ) {
        int edge_id = return_edge_id ( i,j,k,0 );
        map_from_edge_id_to_index_in_E_id.insert (  make_pair ( edge_id, E_id.size() ) );
        E_id.push_back ( edge_id );
      }
    for ( int i=box.x0; i< box.x0+box.dx; i++ ) 
      for ( int j=box.y0; j< box.y0+box.dy-1; j++ )  {
        int edge_id = return_edge_id ( i,j,k,1 );
        map_from_edge_id_to_index_in_E_id.insert (  make_pair ( edge_id, E_id.size() ) );
        E_id.push_back ( edge_id );
      }
  }
  for ( int j=box.y0; j< box.y0+box.dy; j++ )
    for ( int i=box.x0; i< box.x0+box.dx; i++ ) 
      for ( int k = 0; k < numLayers-1; k++ )  {
        int edge_id = return_edge_id ( i,j,k,2 );
        map_from_edge_id_to_index_in_E_id.insert (  make_pair ( edge_id, E_id.size() ) );
        E_id.push_back ( edge_id );
      }
  cout << "Built Grid Box " << endl;
  /* end */
  bnets.clear();
  for ( size_t i=0; i < mini_net_id.size(); i++ ) {
    //  read two terminal pins
    vector<int> terminal_pins ( nets[mini_net_id[i]].terminal_pins.begin(), nets[mini_net_id[i]].terminal_pins.end() );
    int s = terminal_pins[0], t = terminal_pins[1];
    if ( !isInBox (s,box) || !isInBox (t,box) )
      continue;
    //  rip-up those mini nets
    for ( set<int>::iterator it = nets[mini_net_id[i]].edges.begin(); it != nets[mini_net_id[i]].edges.end(); it++ ) {
      -- E[*it].virtual_util;
    }
    //  init the candidate set with current routing
    NetRoute  net_route;
    net_route.edges = nets[mini_net_id[i]].edges;
    // 
    CoarseNet bnet;
    bnet.terminal_pins = nets[mini_net_id[i]].terminal_pins;
    bnet.candidate_routes.push_back ( net_route );
    bnets.push_back ( bnet );
  }
  cout << bnets.size() << endl;
  exit(0);

}

void  CGR::checkCongestedEdgeNumberOfBox ( GridBox &box )
{
  int cnt = 0;
  for ( int k = 0; k < numLayers; k++ )   {
    for ( int j=box.y0; j< box.y0+box.dy; j++ )
      for ( int i=box.x0; i< box.x0+box.dx-1; i++ ) {
        int edge_id = return_edge_id ( i,j,k,0 );
        if ( E[edge_id].cap!=0 && E[edge_id].congestion_ratio > TCR )
          all_the_congested_edges.erase ( edge_id );
      }
    for ( int i=box.x0; i< box.x0+box.dx; i++ ) 
      for ( int j=box.y0; j< box.y0+box.dy-1; j++ )  {
        int edge_id = return_edge_id ( i,j,k,1 );
        if ( E[edge_id].cap!=0 && E[edge_id].congestion_ratio > TCR )
          all_the_congested_edges.erase ( edge_id );
      }
  }
}

//  Added on NOV 11
double  CGR::getNewMetric (void) 
{
  double  metric = 0;
  cout << "TCR = " << TCR << endl;
  for ( vector<GridEdge>::iterator itr = E.begin(); itr != E.end(); itr++ ) {
    if ( (*itr).cap==0 || (*itr).type=="VIA" )
      continue;
    if ( (*itr).congestion_ratio >= TCR )
      metric += ( (*itr).congestion_ratio-TCR );
  }
  return  metric;
}

void  CGR::getTop20 (void )
{
  edge_CR.clear();
  for ( vector<GridEdge>::iterator itr_edge = E.begin(); itr_edge != E.end(); itr_edge++ ) 
    if ( (*itr_edge).cap != 0 && (*itr_edge).type!="VIA"   ) {
      edge_CR.push_back ( make_pair ( (*itr_edge).id, (*itr_edge).congestion_ratio ) );
    }
  
  sort ( edge_CR.begin(), edge_CR.end(), compare_pair_int_double );

  double  y = 0.2;

  for ( size_t i=0; i < edge_CR.size() * y ; i++ ) {
    top20.insert ( edge_CR[i].first );
  }
  cout << "# valid edges = " << edge_CR.size() << endl;

//  TCR = edge_CR[edge_CR.size()*y].second;
  unvisited_top20 = top20;

}

int     CGR::getBoxes3 ( void ) 
{
  boxes.clear();
  for ( size_t i=0; i < top20.size(); i++ ) {
    //  do not waste effort on those blocked but not utilized edges
    if ( E[edge_CR[i].first].util==0 )
      continue;
    bool  visited = false;
    for ( size_t j=0; j < boxes.size(); j++ ) {
      if ( isInBox ( E[edge_CR[i].first].srcPoint, boxes[j] ) && isInBox ( E[edge_CR[i].first].snkPoint, boxes[j] ) ) {
        visited = true;
        break;
      }
    }
    if ( visited==true )
      continue;
    boxes.push_back ( buildBoxFrom3DCongestedEdge ( edge_CR[i].first ) );
  }
  return 1;
}

GridBox   CGR::buildBoxFrom3DCongestedEdge ( int edge_id )
{
  //  printEdge ( edge_id );

  //  via is never congested
  assert ( E[edge_id].type != "VIA" );
  int srcPoint = E[edge_id % edgeCount2D ].srcPoint;
  int snkPoint = E[edge_id % edgeCount2D ].snkPoint;
  GridBox   box;
  box.set_values ( V[srcPoint].x, V[srcPoint].y, V[snkPoint].x-V[srcPoint].x+1, V[snkPoint].y-V[srcPoint].y+1 );
  box.util = E[edge_id].util;
  box.blk = E[edge_id].blk;
  box.cap = E[edge_id].cap;
  box.num_bad_edges = 1;
  GridBox &ref_box = box;

  int itr_num = 0;
  int upper_bound_of_box = 100;
  int lower_bound_of_box = 60;
  int upper_bound_of_bad_edges = 500;

  while ( itr_num < lower_bound_of_box ) {
    int status = expandBox2 ( ref_box );
    ++ itr_num;
  }

  while ( box.num_bad_edges<upper_bound_of_bad_edges && itr_num < upper_bound_of_box && (double)(box.blk+box.util)/(double)box.cap > TCR ) {
    int status = expandBox2 ( ref_box );
    ++ itr_num;
    if ( status == 0 ) {
      cout << "ERROR : self increment bounding box" << endl;
      exit(0);
    }
  }
  //  box.printBox();
  //  cout << box.blk << " " << box.util << " " << box.cap << endl;
  //  cout << "BoxAveCr = " << (double)(box.blk+box.util)/(double)box.cap << endl;
  return  box;
}

pair<double, int> CGR::getBoxAverageCongestionRatio ( GridBox &box ) 
{
  double  sum_cr = 0;
  int cnt = 0;
  for ( int z=0; z< numLayers; z++ ) {
    if ( horizCaps[z] != 0 ) {
      for ( int y=box.y0; y < box.y0+box.dy; y++ )
        for ( int x=box.x0; x < box.x0+box.dx-1; x++ )
          if ( E[return_edge_id ( x,y,z,0 )].cap!=0 ) {
            sum_cr += E[return_edge_id(x,y,z,0)].congestion_ratio;
            ++ cnt;
          }
    }
    if ( vertCaps[z] != 0 ) {
      for ( int x=box.x0; x < box.x0+box.dx; x++ )
        for ( int y=box.y0; y < box.y0+box.dy-1; y++ )
          if ( E[return_edge_id ( x,y,z,1 )].cap!=0 ) {
            sum_cr += E[return_edge_id(x,y,z,1)].congestion_ratio;
            ++ cnt;
          }
    }
  }
  return  make_pair ( sum_cr, cnt );
}

int   CGR::expandBox ( GridBox &box, pair<double,int> box_ave_cr ) 
{
  int lowerleft_xcoordinate = box.x0;
  int lowerleft_ycoordinate = box.y0;
  int width = box.dx;
  int height = box.dy;
  int box_util = box.util;
  int box_cap = box.cap;
  string  L_shape_name;

  vector< pair<string,double> > added_L_congestion;
  map < string, double > L_shape_segment_added_cr;
  map < string, int > L_shape_segment_added_cnt;

  //  
  if ( box.x0-1 >= 0 && box.y0-1 >= 0 ) {
    int added_cr = 0;
    int added_cnt = 0;
    //  lower segment
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        if ( horizCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(i, box.y0-1, z)].ne[1]].congestion_ratio;
          ++ added_cnt;
        }
        if ( vertCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(i, box.y0-1, z)].ne[2]].congestion_ratio;
          ++ added_cnt;
        }
      }
    }
    //  left segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        if ( horizCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(box.x0-1, j, z)].ne[0]].congestion_ratio;
          ++ added_cnt;
        }
        if ( vertCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(box.x0-1, j, z)].ne[3]].congestion_ratio;
          ++ added_cnt;
        }
      }
    }
    L_shape_name = "lower left";
    added_L_congestion.push_back ( make_pair ( L_shape_name,  (double)(box_ave_cr.first+added_cr)/(box_ave_cr.second+added_cnt) ) );
    L_shape_segment_added_cr.insert ( make_pair ( L_shape_name, added_cr ) );
    L_shape_segment_added_cnt.insert ( make_pair ( L_shape_name, added_cnt ) );
  }
  //
  if ( box.x0-1 >= 0 && box.y0+box.dy < fullGridY ) {
    int added_cr = 0;
    int added_cnt = 0;
    //  top segment
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        if ( horizCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(i, box.y0+box.dy, z)].ne[1]].congestion_ratio;
          ++ added_cnt;
        }
        if ( vertCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(i, box.y0+box.dy, z)].ne[3]].congestion_ratio;
          ++ added_cnt;
        }
      }
    }
    //  left segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        if ( horizCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(box.x0-1, j, z)].ne[0]].congestion_ratio;
          ++ added_cnt;
        }
        if ( vertCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(box.x0-1, j, z)].ne[2]].congestion_ratio;
          ++ added_cnt;
        }
      }
    }
    L_shape_name = "top left";
    added_L_congestion.push_back ( make_pair ( L_shape_name,  (double)(box_ave_cr.first+added_cr)/(box_ave_cr.second+added_cnt) ) );
    L_shape_segment_added_cr.insert ( make_pair ( L_shape_name, added_cr ) );
    L_shape_segment_added_cnt.insert ( make_pair ( L_shape_name, added_cnt ) );
  }
  //
  if ( box.x0+box.dx < fullGridX && box.y0-1 >= 0 ) {
    int added_cr = 0;
    int added_cnt = 0;
    //  lower segment
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        if ( horizCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(i, box.y0-1, z)].ne[0]].congestion_ratio;
          ++ added_cnt;
        }
        if ( vertCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(i, box.y0-1, z)].ne[2]].congestion_ratio;
          ++ added_cnt;
        }
      }
    }
    //  right segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        if ( horizCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(box.x0+box.dx, j, z)].ne[1]].congestion_ratio;
          ++ added_cnt;
        }
        if ( vertCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(box.x0+box.dx, j, z)].ne[3]].congestion_ratio;
          ++ added_cnt;
        }
      }
    }
    L_shape_name = "lower right";
    added_L_congestion.push_back ( make_pair ( L_shape_name,  (double)(box_ave_cr.first+added_cr)/(box_ave_cr.second+added_cnt) ) );
    L_shape_segment_added_cr.insert ( make_pair ( L_shape_name, added_cr ) );
    L_shape_segment_added_cnt.insert ( make_pair ( L_shape_name, added_cnt ) );
  }
  //
  if ( box.x0+box.dx < fullGridX && box.y0+box.dy < fullGridY ) {
    int added_cr = 0;
    int added_cnt = 0;
    //  top segment
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        if ( horizCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(i, box.y0+box.dy, z)].ne[0]].congestion_ratio;
          ++ added_cnt;
        }
        if ( vertCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(i, box.y0+box.dy, z)].ne[3]].congestion_ratio;
          ++ added_cnt;
        }
      }
    }
    //  right segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        if ( horizCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(box.x0+box.dx, j, z)].ne[1]].congestion_ratio;
          ++ added_cnt;
        }
        if ( vertCaps[z] != 0 ) {
          added_cr += E[V[xyzToGridNetTerminalId(box.x0+box.dx, j, z)].ne[2]].congestion_ratio;
          ++ added_cnt;
        }
      }
    }
    L_shape_name = "top right";
    added_L_congestion.push_back ( make_pair ( L_shape_name,  (double)(box_ave_cr.first+added_cr)/(box_ave_cr.second+added_cnt) ) );
    L_shape_segment_added_cr.insert ( make_pair ( L_shape_name, added_cr ) );
    L_shape_segment_added_cnt.insert ( make_pair ( L_shape_name, added_cnt ) );
  }

  sort ( added_L_congestion.begin(), added_L_congestion.end(), compare_pair_string_double );
  if ( added_L_congestion[0].first == "lower left" ) {
    box.x0 --; box.y0 --;
    box.dx ++; box.dy ++;
  }
  if ( added_L_congestion[0].first == "lower right" ) {
    box.y0 --;
    box.dx ++; box.dy ++;
  }
  if ( added_L_congestion[0].first == "top left" ) {
    box.x0 --; 
    box.dx ++; box.dy ++;
  }
  if ( added_L_congestion[0].first == "top right" ) {
    box.dx ++; box.dy ++;
  }
  return  1;
}

int   CGR::expandBox2 ( GridBox &box ) 
{
  //
  int edge_id;
  string  selected_L_shape;

  int added_util;
  int added_blk;
  int added_cap;
  int added_num_bad_edges;
  int selected_added_util;
  int selected_added_blk;
  int selected_added_cap;
  int selected_added_num_bad_edges;

  double  box_cr = 100000;  //  init with a large congestion ratio
  //  L0
  if ( box.x0-1 >= 0 && box.y0-1 >= 0 ) {
    added_util = 0;
    added_blk = 0;
    added_cap = 0;
    added_num_bad_edges = 0;
    //  lower segment
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        edge_id = V[xyzToGridNetTerminalId(i,box.y0-1,z)].ne[1];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
        edge_id = V[xyzToGridNetTerminalId(i,box.y0-1,z)].ne[2];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
      }
    }
    //  left segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        edge_id = V[xyzToGridNetTerminalId(box.x0-1,j,z)].ne[0];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
        edge_id = V[xyzToGridNetTerminalId(box.x0-1,j,z)].ne[3];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
      }
    }
    if ( ((double)(box.blk+added_blk+box.util+added_util))/(double)(box.cap+added_cap) < box_cr ) {
      box_cr = ((double)(box.blk+added_blk+box.util+added_util))/(double)(box.cap+added_cap);
      selected_L_shape = "lower left";
      selected_added_util = added_util;
      selected_added_blk = added_blk;
      selected_added_cap = added_cap;
      selected_added_num_bad_edges = added_num_bad_edges;
    }
  }
  //  L1
  if ( box.x0-1 >= 0 && box.y0+box.dy < fullGridY ) {
    added_util = 0;
    added_blk = 0;
    added_cap = 0;
    added_num_bad_edges = 0;
    //  top segment
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        edge_id = V[xyzToGridNetTerminalId(i,box.y0+box.dy,z)].ne[1];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
        edge_id = V[xyzToGridNetTerminalId(i,box.y0+box.dy,z)].ne[3];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
      }
    }
    //  left segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        edge_id = V[xyzToGridNetTerminalId(box.x0-1,j,z)].ne[0];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
        edge_id = V[xyzToGridNetTerminalId(box.x0-1,j,z)].ne[2];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
      }
    }
    if ( ((double)(box.blk+added_blk+box.util+added_util))/(double)(box.cap+added_cap) < box_cr ) {
      box_cr = ((double)(box.blk+added_blk+box.util+added_util))/(double)(box.cap+added_cap);
      selected_L_shape = "top left";
      selected_added_util = added_util;
      selected_added_blk = added_blk;
      selected_added_cap = added_cap;
      selected_added_num_bad_edges = added_num_bad_edges;
    }
  }
  //  L2
  if ( box.x0+box.dx < fullGridX && box.y0-1 >= 0 ) {
    added_util = 0;
    added_blk = 0;
    added_cap = 0;
    added_num_bad_edges = 0;
    //  lower segment
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        edge_id = V[xyzToGridNetTerminalId(i,box.y0-1,z)].ne[0];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
        edge_id = V[xyzToGridNetTerminalId(i,box.y0-1,z)].ne[2];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
      }
    }
    //  right segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        edge_id = V[xyzToGridNetTerminalId(box.x0+box.dx,j,z)].ne[1];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
        edge_id = V[xyzToGridNetTerminalId(box.x0+box.dx,j,z)].ne[3];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
      }
    }
    if ( ((double)(box.blk+added_blk+box.util+added_util))/(double)(box.cap+added_cap) < box_cr ) {
      box_cr = ((double)(box.blk+added_blk+box.util+added_util))/(double)(box.cap+added_cap);
      selected_L_shape = "lower right";
      selected_added_util = added_util;
      selected_added_blk = added_blk;
      selected_added_cap = added_cap;
      selected_added_num_bad_edges = added_num_bad_edges;
    }
  }
  //  L3
  if ( box.x0+box.dx < fullGridX && box.y0+box.dy < fullGridY ) {
    added_util = 0;
    added_blk = 0;
    added_cap = 0;
    added_num_bad_edges = 0;
    //  top segment
    for ( int i=box.x0; i<box.x0+box.dx; i++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        edge_id = V[xyzToGridNetTerminalId(i,box.y0+box.dy,z)].ne[0];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
        edge_id = V[xyzToGridNetTerminalId(i,box.y0+box.dy,z)].ne[3];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
      }
    }
    //  right segment
    for ( int j=box.y0; j<box.y0+box.dy; j++ ) {
      for ( int z=0; z < numLayers; z++ ) {
        edge_id = V[xyzToGridNetTerminalId(box.x0+box.dx,j,z)].ne[1];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
        edge_id = V[xyzToGridNetTerminalId(box.x0+box.dx,j,z)].ne[2];
        assert ( edge_id!=-1 );
        if ( E[edge_id].cap!=0 ) {
          added_util += E[edge_id].util;
          added_blk += E[edge_id].blk;
          added_cap += E[edge_id].cap;
          if ( E[edge_id].congestion_ratio>TCR )
            ++ added_num_bad_edges;
        }
      }
    }
    if ( ((double)(box.blk+added_blk+box.util+added_util))/(double)(box.cap+added_cap) < box_cr ) {
      box_cr = ((double)(box.blk+added_blk+box.util+added_util))/(double)(box.cap+added_cap);
      selected_L_shape = "top right";
      selected_added_util = added_util;
      selected_added_blk = added_blk;
      selected_added_cap = added_cap;
      selected_added_num_bad_edges = added_num_bad_edges;
    }
  }

  if ( selected_L_shape == "lower left" ) {
    box.x0 --; box.y0 --;
    box.dx ++; box.dy ++;
  }
  if ( selected_L_shape == "lower right" ) {
    box.y0 --;
    box.dx ++; box.dy ++;
  }
  if ( selected_L_shape == "top left" ) {
    box.x0 --; 
    box.dx ++; box.dy ++;
  }
  if ( selected_L_shape == "top right" ) {
    box.dx ++; box.dy ++;
  }
  box.util += selected_added_util;
  box.blk += selected_added_blk;
  box.cap += selected_added_cap;
  box.num_bad_edges += selected_added_num_bad_edges;
  return  1;
}

pair<GridBox, int> CGR::getNextBox ( int start )
{
  //  find the first edge which belongs to none of selected boxes
  //  build box for this edge
  int index = start;
  GridBox next_box;
  while ( 1 ) {
    bool  visited = false;
    for ( size_t j=0; j < boxes.size(); j++ ) {
      if ( isInBox ( E[edge_CR[index].first].srcPoint, boxes[j] ) && isInBox ( E[edge_CR[index].first].snkPoint, boxes[j] ) ) {
        visited = true;
        break;
      }
    }
    if ( visited==true ) {
      ++index;
      continue;
    }
    next_box = buildBoxFrom3DCongestedEdge ( edge_CR[index].first );
    boxes.push_back ( next_box );
    //    break;
  }
  return make_pair(next_box, index+1);
}

//  Added on NOV 16
//  add some auxilary functions to monitor the solution within a box
double  CGR::getEdgePenalty ( int edge_id )
{
  double  edge_penalty = 0;
  assert ( mj.size()==bj.size() );
  for ( size_t j=0; j < mj.size(); j++ ) {  //  get the value of edge penalty
    double  penalty_for_single_segment = mj[j] * ( E[edge_id].util+E[edge_id].blk ) / ( E[edge_id].cap+0.0001 ) + bj[j];
    edge_penalty = ( penalty_for_single_segment > edge_penalty ) ? penalty_for_single_segment : edge_penalty;
  } 
  return  edge_penalty;
}

double  CGR::getBoxObjValue ( GridBox &box )
{
  int hut_edge_number = 0;
  int edge_number = 0;
  int edge_id;
  double  box_obj_value = 0;
  int tot_util = 0;
  int tot_blk = 0;
  int tot_cap = 0;
  for ( int z=0; z< numLayers; z++ ) {
    if ( horizCaps[z] != 0 ) {
      for ( int y=box.y0; y < box.y0+box.dy; y++ )
        for ( int x=box.x0; x < box.x0+box.dx-1; x++ ) {
          edge_id = return_edge_id ( x,y,z,0 );
          if ( E[edge_id].cap!=0 && E[edge_id].type!="VIA" ) {
            box_obj_value += getEdgePenalty ( edge_id );
            ++ edge_number;
            tot_util += E[edge_id].util;
            tot_blk += E[edge_id].blk;
            tot_cap += E[edge_id].cap;
            if ( E[edge_id].util>0 )
              ++ hut_edge_number;
          }
        }
    }
    if ( vertCaps[z] != 0 ) {
      for ( int x=box.x0; x < box.x0+box.dx; x++ )
        for ( int y=box.y0; y < box.y0+box.dy-1; y++ ) {
          edge_id = return_edge_id ( x,y,z,1 );
          if ( E[edge_id].cap!=0 && E[edge_id].type!="VIA" ) {
            box_obj_value += getEdgePenalty ( edge_id );
            ++ edge_number;
            tot_util += E[edge_id].util;
            tot_blk += E[edge_id].blk;
            tot_cap += E[edge_id].cap;
            if ( E[edge_id].util>0 )
              ++ hut_edge_number;
          }
        }
    }
  }
  cout << "hut edges " << hut_edge_number << endl;
  cout << "tot edges " << edge_number << endl;
  cout << "total util " << tot_util << endl;
  cout << "total blk " << tot_blk << endl;
  cout << "total cap " << tot_cap << endl;
  return  box_obj_value;
}
//  
//  Added on NOV 17
int     CGR::getNextEdgeToBuildBox (void)
{
  while ( 1 ) {

    bool  visited_in_some_box = false;

    if ( starting_index == edge_CR.size() * 0.2 ) {
      cout << "Done" << endl;
      exit(0);
    }
    for ( size_t j=0; j < boxes.size(); j++ ) {
      if ( isInBox ( E[edge_CR[starting_index].first].srcPoint, boxes[j] ) && isInBox ( E[edge_CR[starting_index].first].snkPoint, boxes[j] ) ) {
        visited_in_some_box = true;
        break;
      }
    }
    if ( visited_in_some_box ==true ) {
      ++ starting_index;
      continue;
    }
    break;
  }
  cout << starting_index << " " << edge_CR.size() * 0.2 << endl;
  return edge_CR[starting_index++].first;
}

void  CGR::buildGridinBox ( GridBox &box )
{
  /* put all the edges inside bounding box into box.E_id */
  E_id.clear();
  box.num_bad_edges = 0;

  congested_edges_number_in_box = 0;
  for ( int k = 0; k < numLayers; k++ )   {
    for ( int j=box.y0; j< box.y0+box.dy; j++ )
      for ( int i=box.x0; i< box.x0+box.dx-1; i++ ) {
        int edge_id = return_edge_id ( i,j,k,0 );
        if (E[edge_id].cap!=0 && E[edge_id].util!=0 && E[edge_id].type!="VIA" && E[edge_id].congestion_ratio>TCR ) {
          all_the_congested_edges.erase(edge_id);
          ++ congested_edges_number_in_box ;
        }
        map_from_edge_id_to_index_in_E_id.insert (  make_pair ( edge_id, E_id.size() ) );
        E_id.push_back ( edge_id );
      }
    for ( int i=box.x0; i< box.x0+box.dx; i++ ) 
      for ( int j=box.y0; j< box.y0+box.dy-1; j++ )  {
        int edge_id = return_edge_id ( i,j,k,1 );
        if (E[edge_id].cap!=0 && E[edge_id].util!=0 && E[edge_id].type!="VIA" && E[edge_id].congestion_ratio>TCR ) {
          all_the_congested_edges.erase(edge_id);
          ++ congested_edges_number_in_box ;
        }
        map_from_edge_id_to_index_in_E_id.insert (  make_pair ( edge_id, E_id.size() ) );
        E_id.push_back ( edge_id );
      }
  }
  for ( int j=box.y0; j< box.y0+box.dy; j++ )
    for ( int i=box.x0; i< box.x0+box.dx; i++ ) 
      for ( int k = 0; k < numLayers-1; k++ )  {
        int edge_id = return_edge_id ( i,j,k,2 );
        map_from_edge_id_to_index_in_E_id.insert (  make_pair ( edge_id, E_id.size() ) );
        E_id.push_back ( edge_id );
      }
  cout << "Built Grid Box " << endl;
  cout << E_id.size() << endl;

  //  compute the threshold in this box after building the grid for box
  vector<double>  edge_congestion_ratio;
  for ( vector<int>::iterator itr_edge = E_id.begin(); itr_edge != E_id.end(); itr_edge++ ) 
    if ( E[*itr_edge].cap!=0 && E[*itr_edge].type!="VIA" ) {
      if ( E[*itr_edge].congestion_ratio > TCR )
        ++ box.num_bad_edges;
      edge_congestion_ratio.push_back ( E[*itr_edge].congestion_ratio );
    }
  sort ( edge_congestion_ratio.begin(), edge_congestion_ratio.end() );
  reverse ( edge_congestion_ratio.begin(), edge_congestion_ratio.end() );

  box.threshold = edge_congestion_ratio[ edge_congestion_ratio.size()*0.2 ];
}

void  CGR::getCongestedNets ( std::ofstream &stream_file )
{
  //  sort edge weights
  vector<pair<int, double> > edge_weight;
  for ( vector<int>::iterator itr = E_id.begin(); itr!= E_id.end(); itr++ )
    if ( E[*itr].type!="VIA" && E[*itr].congestion_ratio>TCR ) {
      edge_weight.push_back ( make_pair( *itr, E[*itr].congestion_ratio ) );
    }
  sort ( edge_weight.begin(), edge_weight.end(), compare_pair_int_double );

  stream_file << "number of bad edges in box = " << edge_weight.size() << endl;;
  stream_file << "TCR = " << TCR << endl;

  considered_nets.clear();

  for ( vector<pair<int,double> >::iterator itr = edge_weight.begin(); itr!= edge_weight.end(); itr++ ) {
    int edge_id = (*itr).first;
    for ( set<int>::iterator itr = E[edge_id].nets_passing_through.begin(); itr != E[edge_id].nets_passing_through.end(); itr++ ) 
      considered_nets.insert ( *itr );
    if ( considered_nets.size() > 200 )
      break;
  }
}

bool  CGR::isInBoxNet ( CoarseNet &net, GridBox &box )
{
  assert ( net.sg.size()!=0 );
  for ( vector<Segment>::iterator it = net.sg.begin(); it != net.sg.end(); it++ )
    if ( !isInBox ( (*it).s, box ) || !isInBox ( (*it).t, box ) ) {
      return  false;
    }
  return  true;
}

void    CGR::genModifiedNetsInBox3 (GridBox &box) 
{
  /* put all the edges inside bounding box into box.E_id */
  E_id.clear();
  bool  debug_mode = false;
  for ( int k = 0; k < numLayers; k++ )   {
    for ( int j=box.y0; j< box.y0+box.dy; j++ )
      for ( int i=box.x0; i< box.x0+box.dx-1; i++ ) {
        int edge_id = return_edge_id ( i,j,k,0 );
        map_from_edge_id_to_index_in_E_id.insert (  make_pair ( edge_id, E_id.size() ) );
        E_id.push_back ( edge_id );
      }
    for ( int i=box.x0; i< box.x0+box.dx; i++ ) 
      for ( int j=box.y0; j< box.y0+box.dy-1; j++ )  {
        int edge_id = return_edge_id ( i,j,k,1 );
        map_from_edge_id_to_index_in_E_id.insert (  make_pair ( edge_id, E_id.size() ) );
        E_id.push_back ( edge_id );
      }
  }
  for ( int j=box.y0; j< box.y0+box.dy; j++ )
    for ( int i=box.x0; i< box.x0+box.dx; i++ ) 
      for ( int k = 0; k < numLayers-1; k++ )  {
        int edge_id = return_edge_id ( i,j,k,2 );
        map_from_edge_id_to_index_in_E_id.insert (  make_pair ( edge_id, E_id.size() ) );
        E_id.push_back ( edge_id );
      }
  cout << "Built Grid Box " << endl;
  cout << E_id.size() << endl;
  /* end */
  bnets.clear();
  for ( set<int>::iterator itr = small_considered_nets.begin(); itr != small_considered_nets.end(); itr++ ) {
    CoarseNet bnet = nets[*itr];
    bnet.old_index = *itr;    //  the old_index in bnet records the original net index
    NetRoute  init_route;
    init_route.edges = bnet.edges;
    bnet.candidate_routes.push_back ( init_route );
    //  remove the corresponding edge util
    for ( set<int>::iterator itr_edge = bnet.edges.begin(); itr_edge != bnet.edges.end(); itr_edge++ ) {
      -- E[*itr_edge].util;
      E[*itr_edge].nets_passing_through.erase ( *itr );
    }
    bnets.push_back ( bnet );
  }
}

void  CGR::readyForNextIteration2 (GridBox &box ) 
{
  for ( vector<CoarseNet>::iterator itr_net = bnets.begin(); itr_net != bnets.end(); itr_net ++ ) {
    //  always keep the candidate routes
    genCandidateRoutingForSimpleNet ( box, *itr_net, 2 );
  }
}

void  CGR::genCandidateRoutingForSimpleNet ( GridBox &box, CoarseNet &net, int M3 )
{
  vector<int> terminal_pins ( net.terminal_pins.begin(), net.terminal_pins.end() );
  int s = terminal_pins[0], t = terminal_pins[1];
  if ( V[s].x==V[t].x || V[s].y==V[t].y )
    genUShapeRoutings ( box, net, s, t, M3 );
  else
    genZShapeRoutings ( box, net, s, t, M3 );
}

void  CGR::genUShapeRoutings ( GridBox &box, CoarseNet &net, int s, int t, int M3 ) 
{
  vector<pair<NetRoute, double> > candidate_routes_weight;
  assert ( s != t );

  int range = 6;
  int c1,c2,c3,c4,c5,c6;
  if ( V[s].y==V[t].y ) {
    for ( int z1 = 0; z1 < numLayers; z1++ ) {
      if ( horizCaps[z1] == 0 )
        continue;
      for ( int y = max(box.y0, V[s].y-range) ; y < min(box.y0+box.dy, V[s].y+range); y++ ) {
        c1 = xyzToGridNetTerminalId ( V[s].x, y, z1 );
        c2 = xyzToGridNetTerminalId ( V[t].x, y, z1 );
        for ( int z2 = 0; z2 < numLayers; z2++ ) {
          if ( vertCaps[z2] == 0 )
            continue;
          NetRoute candidate_route;

          c3 = xyzToGridNetTerminalId ( V[s].x, V[s].y, z2 );
          c4 = xyzToGridNetTerminalId ( V[t].x, V[t].y, z2 );
          c5 = xyzToGridNetTerminalId ( V[s].x, y, z2 );
          c6 = xyzToGridNetTerminalId ( V[t].x, y, z2 );
          set<int>  flat_segment_1 = getFlatSegment ( s, c3 );
          set<int>  flat_segment_2 = getFlatSegment ( c3, c5 );
          set<int>  flat_segment_3 = getFlatSegment ( c5, c1 );
          set<int>  flat_segment_4 = getFlatSegment ( c1, c2 );
          set<int>  flat_segment_5 = getFlatSegment ( c2, c6 );
          set<int>  flat_segment_6 = getFlatSegment ( c6, c4 );
          set<int>  flat_segment_7 = getFlatSegment ( c4, t );
          candidate_route.edges.insert ( flat_segment_1.begin(), flat_segment_1.end() );
          candidate_route.edges.insert ( flat_segment_2.begin(), flat_segment_2.end() );
          candidate_route.edges.insert ( flat_segment_3.begin(), flat_segment_3.end() );
          candidate_route.edges.insert ( flat_segment_4.begin(), flat_segment_4.end() );
          candidate_route.edges.insert ( flat_segment_5.begin(), flat_segment_5.end() );
          candidate_route.edges.insert ( flat_segment_6.begin(), flat_segment_6.end() );
          candidate_route.edges.insert ( flat_segment_7.begin(), flat_segment_7.end() );
          for ( set<int>::iterator it = candidate_route.edges.begin(); it!= candidate_route.edges.end(); it++ )
            assert ( E[*it].cap!=0 );
          candidate_routes_weight.push_back ( make_pair ( candidate_route, computeNetRouteWeight (candidate_route) ) );
        }
      }
    }
  }
  if ( V[s].x==V[t].x  ) {
    for ( int z1 = 0; z1 < numLayers; z1++ ) {
      if ( vertCaps[z1] == 0 )
        continue;
      for ( int x = max(V[s].x-range, box.x0); x < min(V[s].x+range, box.x0+box.dx); x++ ) {
        c1 = xyzToGridNetTerminalId ( x, V[s].y, z1 );
        c2 = xyzToGridNetTerminalId ( x, V[t].y, z1 );
        for ( int z2 = 0; z2 < numLayers; z2++ ) {
          if ( horizCaps[z2] == 0 )
            continue;
          NetRoute  candidate_route;

          c3 = xyzToGridNetTerminalId ( V[s].x, V[s].y, z2 );
          c4 = xyzToGridNetTerminalId ( V[t].x, V[t].y, z2 );
          c5 = xyzToGridNetTerminalId ( x, V[s].y, z2 );
          c6 = xyzToGridNetTerminalId ( x, V[t].y, z2 );
          set<int>  flat_segment_1 = getFlatSegment ( s, c3 );
          set<int>  flat_segment_2 = getFlatSegment ( c3, c5 );
          set<int>  flat_segment_3 = getFlatSegment ( c5, c1 );
          set<int>  flat_segment_4 = getFlatSegment ( c1, c2 );
          set<int>  flat_segment_5 = getFlatSegment ( c2, c6 );
          set<int>  flat_segment_6 = getFlatSegment ( c6, c4 );
          set<int>  flat_segment_7 = getFlatSegment ( c4, t );
          candidate_route.edges.insert ( flat_segment_1.begin(), flat_segment_1.end() );
          candidate_route.edges.insert ( flat_segment_2.begin(), flat_segment_2.end() );
          candidate_route.edges.insert ( flat_segment_3.begin(), flat_segment_3.end() );
          candidate_route.edges.insert ( flat_segment_4.begin(), flat_segment_4.end() );
          candidate_route.edges.insert ( flat_segment_5.begin(), flat_segment_5.end() );
          candidate_route.edges.insert ( flat_segment_6.begin(), flat_segment_6.end() );
          candidate_route.edges.insert ( flat_segment_7.begin(), flat_segment_7.end() );
          for ( set<int>::iterator it = candidate_route.edges.begin(); it!= candidate_route.edges.end(); it++ )
            assert ( E[*it].cap!=0 );
          candidate_routes_weight.push_back ( make_pair ( candidate_route, computeNetRouteWeight (candidate_route) ) );
        }
      }
    }
  }

  sort ( candidate_routes_weight.begin(), candidate_routes_weight.end(), compare_netroute_with_smaller_weight );

  for ( vector<pair<NetRoute, double> >::iterator itr = candidate_routes_weight.begin(); itr < candidate_routes_weight.begin()+M3 && itr != candidate_routes_weight.end(); itr++ )
    net.candidate_routes.push_back ( (*itr).first );
}

void  CGR::genZShapeRoutings ( GridBox &box, CoarseNet &net, int s, int t, int M3 ) 
{
  vector<pair<NetRoute, double> > candidate_routes_weight;
  assert ( s != t );

  int c1,c2,c3,c4,c5,c6;
  //  
  for ( int y = min(V[s].y, V[t].y); y <= max(V[s].y, V[t].y); y++ ) {
    for ( int z1=0; z1< numLayers; z1++ ) {
      if ( vertCaps[z1]==0 )
        continue;
      for ( int z2=0; z2< numLayers; z2++ ) {
        if ( horizCaps[z2]==0 )
          continue;
        for ( int z3=0; z3< numLayers; z3++ ) {
          if ( vertCaps[z3]==0 )
            continue;
          NetRoute  candidate_route;

          c1 = xyzToGridNetTerminalId ( V[s].x, V[s].y, z1 );
          c2 = xyzToGridNetTerminalId ( V[s].x, y, z1 );
          c3 = xyzToGridNetTerminalId ( V[s].x, y, z2 );
          c4 = xyzToGridNetTerminalId ( V[t].x, y, z2 );
          c5 = xyzToGridNetTerminalId ( V[t].x, y, z3 );
          c6 = xyzToGridNetTerminalId ( V[t].x, V[t].y, z3 );

          set<int>  flat_segment_1 = getFlatSegment ( s, c1 );
          set<int>  flat_segment_2 = getFlatSegment ( c1, c2 );
          set<int>  flat_segment_3 = getFlatSegment ( c2, c3 );
          set<int>  flat_segment_4 = getFlatSegment ( c3, c4 );
          set<int>  flat_segment_5 = getFlatSegment ( c4, c5 );
          set<int>  flat_segment_6 = getFlatSegment ( c5, c6 );
          set<int>  flat_segment_7 = getFlatSegment ( c6, t );
          candidate_route.edges.insert ( flat_segment_1.begin(), flat_segment_1.end() );
          candidate_route.edges.insert ( flat_segment_2.begin(), flat_segment_2.end() );
          candidate_route.edges.insert ( flat_segment_3.begin(), flat_segment_3.end() );
          candidate_route.edges.insert ( flat_segment_4.begin(), flat_segment_4.end() );
          candidate_route.edges.insert ( flat_segment_5.begin(), flat_segment_5.end() );
          candidate_route.edges.insert ( flat_segment_6.begin(), flat_segment_6.end() );
          candidate_route.edges.insert ( flat_segment_7.begin(), flat_segment_7.end() );
          for ( set<int>::iterator it = candidate_route.edges.begin(); it!= candidate_route.edges.end(); it++ )
            assert ( E[*it].cap!=0 );
          candidate_routes_weight.push_back ( make_pair ( candidate_route, computeNetRouteWeight (candidate_route) ) );
        }
      }
    }
  }
  //
  for ( int x = min(V[s].x, V[t].x); x <= max(V[s].x, V[t].x); x++ ) {
    for ( int z1=0; z1< numLayers; z1++ ) {
      if ( horizCaps[z1]==0 )
        continue;
      for ( int z2=0; z2< numLayers; z2++ ) {
        if ( vertCaps[z2]==0 )
          continue;
        for ( int z3=0; z3< numLayers; z3++ ) {
          if ( horizCaps[z3]==0 )
            continue;
          NetRoute  candidate_route;

          c1 = xyzToGridNetTerminalId ( V[s].x, V[s].y, z1 );
          c2 = xyzToGridNetTerminalId ( x, V[s].y, z1 );
          c3 = xyzToGridNetTerminalId ( x, V[s].y, z2 );
          c4 = xyzToGridNetTerminalId ( x, V[t].y, z2 );
          c5 = xyzToGridNetTerminalId ( x, V[t].y, z3 );
          c6 = xyzToGridNetTerminalId ( V[t].x, V[t].y, z3 );

          set<int>  flat_segment_1 = getFlatSegment ( s, c1 );
          set<int>  flat_segment_2 = getFlatSegment ( c1, c2 );
          set<int>  flat_segment_3 = getFlatSegment ( c2, c3 );
          set<int>  flat_segment_4 = getFlatSegment ( c3, c4 );
          set<int>  flat_segment_5 = getFlatSegment ( c4, c5 );
          set<int>  flat_segment_6 = getFlatSegment ( c5, c6 );
          set<int>  flat_segment_7 = getFlatSegment ( c6, t );
          candidate_route.edges.insert ( flat_segment_1.begin(), flat_segment_1.end() );
          candidate_route.edges.insert ( flat_segment_2.begin(), flat_segment_2.end() );
          candidate_route.edges.insert ( flat_segment_3.begin(), flat_segment_3.end() );
          candidate_route.edges.insert ( flat_segment_4.begin(), flat_segment_4.end() );
          candidate_route.edges.insert ( flat_segment_5.begin(), flat_segment_5.end() );
          candidate_route.edges.insert ( flat_segment_6.begin(), flat_segment_6.end() );
          candidate_route.edges.insert ( flat_segment_7.begin(), flat_segment_7.end() );
          for ( set<int>::iterator it = candidate_route.edges.begin(); it!= candidate_route.edges.end(); it++ )
            assert ( E[*it].cap!=0 );
          candidate_routes_weight.push_back ( make_pair ( candidate_route, computeNetRouteWeight (candidate_route) ) );
        }
      }
    }
  }
  /* Select the top M2 candidate subnets with least weight */
  if ( M3==2 ) {
    if ( candidate_routes_weight.size()==0 )
      return;
    vector<pair<NetRoute, double> >::iterator itr_smallest = candidate_routes_weight.begin();
    for ( vector<pair<NetRoute, double> >::iterator it = candidate_routes_weight.begin(); it != candidate_routes_weight.end(); it++ ) 
      if ( (*it).second < (*itr_smallest).second ) 
        itr_smallest = it ;
    net.candidate_routes.push_back ( (*itr_smallest).first );

    return;
  }
  sort ( candidate_routes_weight.begin(), candidate_routes_weight.end(), compare_netroute_with_smaller_weight );

  for ( vector<pair<NetRoute, double> >::iterator itr = candidate_routes_weight.begin(); itr < candidate_routes_weight.begin()+M3 && itr != candidate_routes_weight.end(); itr++ )
    net.candidate_routes.push_back ( (*itr).first );
}

//
//  Added on NOV 19
//  box selection algorithm
void  CGR::selectBox (void)
{
  vector<int> hut_edges_to_fix ( involved_hut_edges.begin(), involved_hut_edges.end() );

  int current_index = 0;
  while ( current_index < hut_edges_to_fix.size() ) {
    bool  visited_in_some_box = false;
    for ( size_t j=0; j < boxes.size(); j++ ) {
      if ( isInBox ( E[hut_edges_to_fix[current_index]].srcPoint, boxes[j] ) && isInBox ( E[hut_edges_to_fix[current_index]].snkPoint, boxes[j] ) ) {
        visited_in_some_box = true;
        break;    //  break the for loop
      }
    }
    if ( visited_in_some_box == false) {
      GridBox b = buildBoxFrom3DCongestedEdge ( hut_edges_to_fix[current_index] );
      boxes.push_back ( b );
    }
    ++ current_index ;
  }
}

void    CGR::genModifiedNetsInBox4 (GridBox &box) 
{
  /* put all the edges inside bounding box into box.E_id */
  E_id.clear();
  bool  debug_mode = true;
  for ( int k = 0; k < numLayers; k++ )   {
    for ( int j=box.y0; j< box.y0+box.dy; j++ )
      for ( int i=box.x0; i< box.x0+box.dx-1; i++ ) {
        int edge_id = return_edge_id ( i,j,k,0 );
        map_from_edge_id_to_index_in_E_id.insert (  make_pair ( edge_id, E_id.size() ) );
        E_id.push_back ( edge_id );
      }
    for ( int i=box.x0; i< box.x0+box.dx; i++ ) 
      for ( int j=box.y0; j< box.y0+box.dy-1; j++ )  {
        int edge_id = return_edge_id ( i,j,k,1 );
        map_from_edge_id_to_index_in_E_id.insert (  make_pair ( edge_id, E_id.size() ) );
        E_id.push_back ( edge_id );
      }
  }
  for ( int j=box.y0; j< box.y0+box.dy; j++ )
    for ( int i=box.x0; i< box.x0+box.dx; i++ ) 
      for ( int k = 0; k < numLayers-1; k++ )  {
        int edge_id = return_edge_id ( i,j,k,2 );
        map_from_edge_id_to_index_in_E_id.insert (  make_pair ( edge_id, E_id.size() ) );
        E_id.push_back ( edge_id );
      }
  cout << "Built Grid Box " << endl;
  cout << E_id.size() << endl;
  /* end */
  bnets.clear();
  for ( set<int>::iterator itr = considered_nets.begin(); itr != considered_nets.end(); itr++ ) {
    CoarseNet bnet = genSingleModifiedNetInBox ( nets[*itr], box );
    bnet.old_index = *itr;    //  the old_index in bnet records the original net index
    NetRoute  init_route;
    init_route.edges = bnet.edges;
    bnet.candidate_routes.push_back ( init_route );
    //  remove the corresponding edge util
    for ( set<int>::iterator itr_edge = bnet.edges.begin(); itr_edge != bnet.edges.end(); itr_edge++ ) {
      -- E[*itr_edge].util;
      E[*itr_edge].nets_passing_through.erase ( *itr );
    }
    netDecomposition ( bnet );
    if ( debug_mode == true ) {
      if ( isConnected ( bnet ) == false )
        cout << bnet.old_index << " disconnected" << endl;
      else {
        cout << bnet.old_index << " " << bnet.terminal_pins.size() << " " << bnet.edges.size() << " " << bnet.sg.size() << " " << bnet.subnets.size() << endl;
        for ( size_t ii=0; ii < bnet.subnets.size(); ii++ ) {
          cout << "\t" << bnet.subnets[ii].edges.size() <<"(";
          bnet.subnets[ii].is_congested = false;
          for ( set<int>::iterator itr = bnet.subnets[ii].edges.begin(); itr != bnet.subnets[ii].edges.end(); itr++ )
            if ( box.congested_edges.find(*itr)!= box.congested_edges.end() ) {
              bnet.subnets[ii].is_congested = true;
              break;
            }
          if ( bnet.subnets[ii].is_congested == true )
            cout << "Y)";
          else
            cout << "N)";
        }
        cout << endl;
      }
    }
    bnets.push_back ( bnet );
  }
  cout << bnets.size() << endl;
}

//  Added on NOV 21
void  CGR::readyForNextIteration3 (GridBox &box ) 
{
  for ( size_t i=0; i< bnets.size(); i++ )
    genCandidateRoute2 ( box, i, 2 );
}

void  CGR::genCandidateRoute (GridBox &box, int bnet_index, int M3 )
{
  CoarseNet &bnet = bnets[bnet_index];
  if ( isConnected(bnet)==false )
    return;
  if ( bnet.subnets.size()==1 && getNetHPWL ( bnet ) > 10 )
    return;

  //  initialize candidate subnets
  if ( bnet.candidate_subnets.size()==0 ) {
    for ( vector<TwoTerCoarseNet>::iterator itr_subnet = bnet.subnets.begin(); itr_subnet != bnet.subnets.end(); ++itr_subnet ) {
      vector<TwoTerCoarseNet> snet_routes;
      snet_routes.push_back ( *itr_subnet );
      bnet.candidate_subnets.push_back ( snet_routes );
    }
    assert ( bnet.subnets.size()==bnet.candidate_subnets.size() );
  }
  else 
    for ( size_t j=0; j<bnet.subnets.size(); j++ ) {
      bnet.candidate_subnets[j].clear();
      bnet.candidate_subnets[j].push_back( bnet.subnets[j] );
    }

  //  generate candidate routes for subnets
  for ( size_t j=0; j<bnet.subnets.size(); j++ ) {
    TwoTerCoarseNet &snet = bnet.subnets[j];
    if ( bnet.subnets[j].is_congested == false )
      continue;
    //  generate candidate routes for subnet
    genCandidateRoutingForSimpleSubnet ( box, make_pair(bnet_index,j), M3 );
  }
  addNewCandidateRoutes ( bnet_index, M3 );

  bool  debug_mode = false;
  if ( debug_mode == true ) {
    //  print the results
    cout << bnet.old_index << " " << bnet.terminal_pins.size() << " " << bnet.edges.size() << " " << bnet.sg.size() << " " << bnet.subnets.size() << " " << getNetHPWL(bnet) << endl;
    for ( size_t ii=0; ii < bnet.subnets.size(); ii++ ) {
      cout << "\t" << bnet.subnets[ii].edges.size() <<"(";
      bnet.subnets[ii].is_congested = false;
      for ( set<int>::iterator itr = bnet.subnets[ii].edges.begin(); itr != bnet.subnets[ii].edges.end(); itr++ )
        if ( box.congested_edges.find(*itr)!= box.congested_edges.end() ) {
          bnet.subnets[ii].is_congested = true;
          break;
        }
      if ( bnet.subnets[ii].is_congested == true )
        cout << "Y " << bnet.candidate_subnets[ii].size() <<") ";
      else
        cout << "N " << bnet.candidate_subnets[ii].size() <<") ";
    }
    cout << endl;
  }

}

void  CGR::genCandidateRoutingForSimpleSubnet ( GridBox &box, TwoTerCoarseNetIdType subnet_id , int M3 )
{
  int ii = subnet_id.first;
  int jj = subnet_id.second;
  int s = bnets[ii].subnets[jj].s;
  int t = bnets[ii].subnets[jj].t;
  if ( V[s].x==V[t].x || V[s].y==V[t].y )
    genUShapeRoutingsForSubnet ( box, subnet_id, M3 );
  else
    genZShapeRoutingsForSubnet ( box, subnet_id, M3 );
}

void  CGR::genUShapeRoutingsForSubnet ( GridBox &box, TwoTerCoarseNetIdType subnet_id, int M3 )
{
  vector<pair<TwoTerCoarseNet, double> > candidate_routes_weight;
  int ii = subnet_id.first;
  int jj = subnet_id.second;
  int s = bnets[ii].subnets[jj].s;
  int t = bnets[ii].subnets[jj].t;
  assert ( s != t );

  int range = 6;
  int c1,c2,c3,c4,c5,c6;
  if ( V[s].y==V[t].y ) {
    for ( int z1 = 0; z1 < numLayers; z1++ ) {
      if ( horizCaps[z1] == 0 )
        continue;
      for ( int y = max(box.y0, V[s].y-range) ; y < min(box.y0+box.dy, V[s].y+range); y++ ) {
        c1 = xyzToGridNetTerminalId ( V[s].x, y, z1 );
        c2 = xyzToGridNetTerminalId ( V[t].x, y, z1 );
        for ( int z2 = 0; z2 < numLayers; z2++ ) {
          if ( vertCaps[z2] == 0 )
            continue;
          TwoTerCoarseNet candidate_route;
          candidate_route.s = s;
          candidate_route.t = t;

          c3 = xyzToGridNetTerminalId ( V[s].x, V[s].y, z2 );
          c4 = xyzToGridNetTerminalId ( V[t].x, V[t].y, z2 );
          c5 = xyzToGridNetTerminalId ( V[s].x, y, z2 );
          c6 = xyzToGridNetTerminalId ( V[t].x, y, z2 );
          set<int>  flat_segment_1 = getFlatSegment ( s, c3 );
          set<int>  flat_segment_2 = getFlatSegment ( c3, c5 );
          set<int>  flat_segment_3 = getFlatSegment ( c5, c1 );
          set<int>  flat_segment_4 = getFlatSegment ( c1, c2 );
          set<int>  flat_segment_5 = getFlatSegment ( c2, c6 );
          set<int>  flat_segment_6 = getFlatSegment ( c6, c4 );
          set<int>  flat_segment_7 = getFlatSegment ( c4, t );
          candidate_route.edges.insert ( flat_segment_1.begin(), flat_segment_1.end() );
          candidate_route.edges.insert ( flat_segment_2.begin(), flat_segment_2.end() );
          candidate_route.edges.insert ( flat_segment_3.begin(), flat_segment_3.end() );
          candidate_route.edges.insert ( flat_segment_4.begin(), flat_segment_4.end() );
          candidate_route.edges.insert ( flat_segment_5.begin(), flat_segment_5.end() );
          candidate_route.edges.insert ( flat_segment_6.begin(), flat_segment_6.end() );
          candidate_route.edges.insert ( flat_segment_7.begin(), flat_segment_7.end() );
          double  route_weight = 0;
          for ( set<int>::iterator it = candidate_route.edges.begin(); it!= candidate_route.edges.end(); it++ )
            route_weight += E[*it].congestion_ratio;
          candidate_routes_weight.push_back ( make_pair ( candidate_route, route_weight ) );
        }
      }
    }
  }
  if ( V[s].x==V[t].x  ) {
    for ( int z1 = 0; z1 < numLayers; z1++ ) {
      if ( vertCaps[z1] == 0 )
        continue;
      for ( int x = max(V[s].x-range, box.x0); x < min(V[s].x+range, box.x0+box.dx); x++ ) {
        c1 = xyzToGridNetTerminalId ( x, V[s].y, z1 );
        c2 = xyzToGridNetTerminalId ( x, V[t].y, z1 );
        for ( int z2 = 0; z2 < numLayers; z2++ ) {
          if ( horizCaps[z2] == 0 )
            continue;
          TwoTerCoarseNet candidate_route;
          candidate_route.s = s;
          candidate_route.t = t;

          c3 = xyzToGridNetTerminalId ( V[s].x, V[s].y, z2 );
          c4 = xyzToGridNetTerminalId ( V[t].x, V[t].y, z2 );
          c5 = xyzToGridNetTerminalId ( x, V[s].y, z2 );
          c6 = xyzToGridNetTerminalId ( x, V[t].y, z2 );
          set<int>  flat_segment_1 = getFlatSegment ( s, c3 );
          set<int>  flat_segment_2 = getFlatSegment ( c3, c5 );
          set<int>  flat_segment_3 = getFlatSegment ( c5, c1 );
          set<int>  flat_segment_4 = getFlatSegment ( c1, c2 );
          set<int>  flat_segment_5 = getFlatSegment ( c2, c6 );
          set<int>  flat_segment_6 = getFlatSegment ( c6, c4 );
          set<int>  flat_segment_7 = getFlatSegment ( c4, t );
          candidate_route.edges.insert ( flat_segment_1.begin(), flat_segment_1.end() );
          candidate_route.edges.insert ( flat_segment_2.begin(), flat_segment_2.end() );
          candidate_route.edges.insert ( flat_segment_3.begin(), flat_segment_3.end() );
          candidate_route.edges.insert ( flat_segment_4.begin(), flat_segment_4.end() );
          candidate_route.edges.insert ( flat_segment_5.begin(), flat_segment_5.end() );
          candidate_route.edges.insert ( flat_segment_6.begin(), flat_segment_6.end() );
          candidate_route.edges.insert ( flat_segment_7.begin(), flat_segment_7.end() );
          double  route_weight = 0;
          for ( set<int>::iterator it = candidate_route.edges.begin(); it!= candidate_route.edges.end(); it++ )
            route_weight += E[*it].congestion_ratio;
          candidate_routes_weight.push_back ( make_pair ( candidate_route, route_weight ) );
        }
      }
    }
  }

  sort ( candidate_routes_weight.begin(), candidate_routes_weight.end(), compare_subnet_with_smaller_weight );

  for ( vector<pair<TwoTerCoarseNet, double> >::iterator itr = candidate_routes_weight.begin(); itr < candidate_routes_weight.begin()+M3 && itr != candidate_routes_weight.end(); itr++ )
    if ( (*itr).first != bnets[ii].subnets[jj] )
      bnets[ii].candidate_subnets[jj].push_back ( (*itr).first );
}

void  CGR::genZShapeRoutingsForSubnet ( GridBox &box, TwoTerCoarseNetIdType subnet_id, int M3 )
{
  vector<pair<TwoTerCoarseNet, double> > candidate_routes_weight;
  int ii = subnet_id.first;
  int jj = subnet_id.second;
  int s = bnets[ii].subnets[jj].s;
  int t = bnets[ii].subnets[jj].t;
  assert ( s != t );

  int c1,c2,c3,c4,c5,c6;
  //  
  for ( int y = min(V[s].y, V[t].y); y <= max(V[s].y, V[t].y); y++ ) {
    for ( int z1=0; z1< numLayers; z1++ ) {
      if ( vertCaps[z1]==0 )
        continue;
      for ( int z2=0; z2< numLayers; z2++ ) {
        if ( horizCaps[z2]==0 )
          continue;
        for ( int z3=0; z3< numLayers; z3++ ) {
          if ( vertCaps[z3]==0 )
            continue;
          TwoTerCoarseNet candidate_route;
          candidate_route.s = s;
          candidate_route.t = t;

          c1 = xyzToGridNetTerminalId ( V[s].x, V[s].y, z1 );
          c2 = xyzToGridNetTerminalId ( V[s].x, y, z1 );
          c3 = xyzToGridNetTerminalId ( V[s].x, y, z2 );
          c4 = xyzToGridNetTerminalId ( V[t].x, y, z2 );
          c5 = xyzToGridNetTerminalId ( V[t].x, y, z3 );
          c6 = xyzToGridNetTerminalId ( V[t].x, V[t].y, z3 );

          set<int>  flat_segment_1 = getFlatSegment ( s, c1 );
          set<int>  flat_segment_2 = getFlatSegment ( c1, c2 );
          set<int>  flat_segment_3 = getFlatSegment ( c2, c3 );
          set<int>  flat_segment_4 = getFlatSegment ( c3, c4 );
          set<int>  flat_segment_5 = getFlatSegment ( c4, c5 );
          set<int>  flat_segment_6 = getFlatSegment ( c5, c6 );
          set<int>  flat_segment_7 = getFlatSegment ( c6, t );
          candidate_route.edges.insert ( flat_segment_1.begin(), flat_segment_1.end() );
          candidate_route.edges.insert ( flat_segment_2.begin(), flat_segment_2.end() );
          candidate_route.edges.insert ( flat_segment_3.begin(), flat_segment_3.end() );
          candidate_route.edges.insert ( flat_segment_4.begin(), flat_segment_4.end() );
          candidate_route.edges.insert ( flat_segment_5.begin(), flat_segment_5.end() );
          candidate_route.edges.insert ( flat_segment_6.begin(), flat_segment_6.end() );
          candidate_route.edges.insert ( flat_segment_7.begin(), flat_segment_7.end() );
          double  route_weight = 0;
          for ( set<int>::iterator it = candidate_route.edges.begin(); it!= candidate_route.edges.end(); it++ )
            route_weight += E[*it].congestion_ratio;
          candidate_routes_weight.push_back ( make_pair ( candidate_route, route_weight ) );
        }
      }
    }
  }
  //
  for ( int x = min(V[s].x, V[t].x); x <= max(V[s].x, V[t].x); x++ ) {
    for ( int z1=0; z1< numLayers; z1++ ) {
      if ( horizCaps[z1]==0 )
        continue;
      for ( int z2=0; z2< numLayers; z2++ ) {
        if ( vertCaps[z2]==0 )
          continue;
        for ( int z3=0; z3< numLayers; z3++ ) {
          if ( horizCaps[z3]==0 )
            continue;
          TwoTerCoarseNet candidate_route;
          candidate_route.s = s;
          candidate_route.t = t;

          c1 = xyzToGridNetTerminalId ( V[s].x, V[s].y, z1 );
          c2 = xyzToGridNetTerminalId ( x, V[s].y, z1 );
          c3 = xyzToGridNetTerminalId ( x, V[s].y, z2 );
          c4 = xyzToGridNetTerminalId ( x, V[t].y, z2 );
          c5 = xyzToGridNetTerminalId ( x, V[t].y, z3 );
          c6 = xyzToGridNetTerminalId ( V[t].x, V[t].y, z3 );

          set<int>  flat_segment_1 = getFlatSegment ( s, c1 );
          set<int>  flat_segment_2 = getFlatSegment ( c1, c2 );
          set<int>  flat_segment_3 = getFlatSegment ( c2, c3 );
          set<int>  flat_segment_4 = getFlatSegment ( c3, c4 );
          set<int>  flat_segment_5 = getFlatSegment ( c4, c5 );
          set<int>  flat_segment_6 = getFlatSegment ( c5, c6 );
          set<int>  flat_segment_7 = getFlatSegment ( c6, t );
          candidate_route.edges.insert ( flat_segment_1.begin(), flat_segment_1.end() );
          candidate_route.edges.insert ( flat_segment_2.begin(), flat_segment_2.end() );
          candidate_route.edges.insert ( flat_segment_3.begin(), flat_segment_3.end() );
          candidate_route.edges.insert ( flat_segment_4.begin(), flat_segment_4.end() );
          candidate_route.edges.insert ( flat_segment_5.begin(), flat_segment_5.end() );
          candidate_route.edges.insert ( flat_segment_6.begin(), flat_segment_6.end() );
          candidate_route.edges.insert ( flat_segment_7.begin(), flat_segment_7.end() );
          double  route_weight = 0;
          for ( set<int>::iterator it = candidate_route.edges.begin(); it!= candidate_route.edges.end(); it++ )
            route_weight += E[*it].congestion_ratio;
          candidate_routes_weight.push_back ( make_pair ( candidate_route, route_weight ) );
        }
      }
    }
  }
  sort ( candidate_routes_weight.begin(), candidate_routes_weight.end(), compare_subnet_with_smaller_weight );

  for ( vector<pair<TwoTerCoarseNet, double> >::iterator itr = candidate_routes_weight.begin(); itr < candidate_routes_weight.begin()+M3 && itr != candidate_routes_weight.end(); itr++ )
    if ( (*itr).first != bnets[ii].subnets[jj] )
      bnets[ii].candidate_subnets[jj].push_back ( (*itr).first );
}

void  CGR::genCandidateRoutingForComplexSubnet ( GridBox &box, TwoTerCoarseNetIdType subnet_id , int M3 )
{
  int ii = subnet_id.first;
  int jj = subnet_id.second;
  int s = bnets[ii].subnets[jj].s;
  int t = bnets[ii].subnets[jj].t;

  if ( V[s].x==V[t].x || V[s].y==V[t].y )
    genUShapeRoutingsForSubnet ( box, subnet_id, M3 );
  else
    genZShapeRoutingsForSubnet ( box, subnet_id, M3 );
}

void    CGR::genModifiedNetsInBox5 (GridBox &box) 
{
  bool  debug_mode = false;
  for ( set<int>::iterator itr = considered_nets.begin(); itr != considered_nets.end(); itr++ ) {
    CoarseNet &net = nets[*itr];
    CoarseNet bnet = genSingleModifiedNetInBox ( net, box );

    assert ( bnet.edges.empty() );
    for ( size_t  i=0;  i< bnet.sg.size(); i++ ) {
      vector<int> edge_vector = getEdgeVectorFromSegment ( bnet.sg[i] );
      for ( size_t j=0; j< edge_vector.size(); j++ ) {
        assert ( net.edges.find(edge_vector[j])!=net.edges.end() );
        bnet.edges.insert ( edge_vector[j] );
      }
    }
    if ( isConnected ( bnet )==false )
      continue;
    for ( set<int>::iterator it = bnet.edges.begin(); it != bnet.edges.end(); it++ ) {
      assert ( net.edges.find(*it)!=net.edges.end() );
      net.edges.erase(*it);
      -- E[*it].util;
      assert ( E[*it].nets_passing_through.find(*itr)!=E[*it].nets_passing_through.end() );
      E[*it].nets_passing_through.erase(*itr);
    }

    bnet.old_index = *itr;    //  the old_index in bnet records the original net index

    netDecomposition ( bnet );

    NetRoute  init_route;
    init_route.edges = bnet.edges;
    bnet.candidate_routes.push_back ( init_route );

    bnets.push_back ( bnet );
  }
}

void  CGR::refreshNetSegments (void) 
{
  for ( size_t i=0; i< bnets.size(); i++ ) {

    int net_id =  bnets[i].old_index;

    for ( set<int>::iterator itr = bnets[i].edges.begin(); itr!= bnets[i].edges.end(); itr++ ) {
      assert ( E[*itr].cap!=0 );
      /* insert those edges inside bounding box using the solution of ILP */
      assert ( nets[net_id].edges.find(*itr)==nets[net_id].edges.end() );
      nets[net_id].edges.insert ( *itr );
      //  update edge utilization
      assert ( E[*itr].nets_passing_through.find(net_id)==E[*itr].nets_passing_through.end() );
      E[*itr].nets_passing_through.insert ( net_id );
      ++ E[*itr].util;
    }
    /* update segment vector in net[index_in_nets] */
    generateLongestSegments ( net_id );
  }

  for ( vector<int>::iterator it = E_id.begin(); it!= E_id.end(); it++ )
    E[*it].congestion_ratio = (double)(E[*it].blk + E[*it].util)/(E[*it].cap+0.0001);
}

void  CGR::readyForNextIteration4 (GridBox &box ) 
{
  //  this function only adds new candidate routes, doesn't change edge util or something else
  cout << "bnets.size()=" << bnets.size() << endl;
  for ( size_t i=0; i< bnets.size(); i++ )  {
    CoarseNet &bnet = bnets[i];
    if ( bnet.candidate_routes.size()==1 && bnet.terminal_pins.size()==2 && getNetHPWL(bnet)<40 )
      genCandidateRoutingForSimpleNet ( box, bnet, 1 );
    if ( bnet.candidate_routes.size()!=1 ) 
      genCandidateRoute4( box, i, 2 );
  }
}

int   CGR::getCongestedEdgeNumberOfBox ( GridBox &box )
{
  //  don't need it
  return  0;
}

void  CGR::genCandidateRoute2 ( GridBox &box, int bnet_index, int M3 )
{
  CoarseNet &bnet = bnets[bnet_index];
  bnet.congested_edges.clear();

  //  identify all the congested edges in current route
  cout << "#" << endl;
  cout << "terminal pins" << endl;
  for ( set<int>::iterator it = bnet.terminal_pins.begin(); it!=bnet.terminal_pins.end(); it++ )
    printNode (*it);
  cout << "edges" << endl;
  for ( set<int>::iterator it = bnet.edges.begin(); it!=bnet.edges.end(); it++ )
    printEdge(*it);
  for ( set<int>::iterator itr = bnet.edges.begin(); itr != bnet.edges.end(); itr++ ) {
    if ( E[*itr].congestion_ratio > TCR ) {
      bnet.congested_edges.insert ( *itr );
    }
  }
  //  reroute each congested edge
  for ( set<int>::iterator itr = bnet.congested_edges.begin(); itr!= bnet.congested_edges.end(); itr++ ) {
    detourSingleEdge ( box, bnet_index, *itr );
  }
  NetRoute new_candidate_route;
  new_candidate_route.edges = bnet.edges;
  bnet.candidate_routes.push_back ( new_candidate_route );
  //  print result
  string  str = bnet.candidate_routes[0].edges!=bnet.candidate_routes[1].edges?"TRUE":"FALSE" ;
  cout << bnet.old_index << " " << bnet.terminal_pins.size() << " " << getNetHPWL ( bnet ) << " " << bnet.congested_edges.size() << " " << bnet.candidate_routes.size() << " " << str << " " << bnet.candidate_routes[0].edges.size() << " " << bnet.candidate_routes[1].edges.size() << endl;
  cout << "route 0" << endl;
  for ( set<int>::iterator it = bnet.candidate_routes[0].edges.begin(); it!=bnet.candidate_routes[0].edges.end(); it++ )
    printEdge(*it);
  cout << "route 1" << endl;
  for ( set<int>::iterator it = bnet.candidate_routes[1].edges.begin(); it!=bnet.candidate_routes[1].edges.end(); it++ )
    printEdge(*it);
}

bool  CGR::isTerminalPin ( CoarseNet &bnet, int node )
{
  return  ( bnet.terminal_pins.find(node)!=bnet.terminal_pins.end() );
}

bool  CGR::isDanglingNode ( CoarseNet &bnet, int node, set<int> &reroute )
{
  if ( isTerminalPin ( bnet, node ) )
    return  false;
  int cnt=0;
  for ( int i=0; i < 5; i++ )
    if ( reroute.find( V[node].ne[i] )!=reroute.end() )
      ++ cnt;
  return  ( cnt==1 );
}

bool  CGR::detourPath ( CoarseNet &bnet, vector<int> node_path, double &detour_cost, set<int> &reroute )
{
  vector<int> edge_path;
  for ( size_t i=0; i < node_path.size()-1; i++ ) {
    int edge_id = getEdgeIdWithTwoEndingPoints ( node_path[i], node_path[i+1] );
    //  debug
    if ( edge_id==-1 ) {
      printNode(node_path[i]);
      printNode(node_path[i+1]);
      cout<< "1" << endl;
      for ( size_t i=0; i< node_path.size(); i++ )
        printNode ( node_path[i] );
      return  false;
    }
    //    printEdge ( edge_id );
    if ( bnet.edges.find(edge_id)==bnet.edges.end() ) {
      reroute.insert ( edge_id );
      detour_cost += E[edge_id].weight;
      //      printEdge ( edge_id );
      //      cout << "+" << E[edge_id].weight << endl;
    }
  }
  int  i=0;
  //  cout << "after 1 " << detour_cost << endl;
  while ( i+1 < node_path.size() && isDanglingNode(bnet, node_path[i], reroute )) {
    int edge_id = getEdgeIdWithTwoEndingPoints ( node_path[i], node_path[i+1] );
    if ( edge_id==-1 ) {
      printNode(node_path[i]);
      printNode(node_path[i+1]);
      cout<< "2" << endl;
      for ( size_t i=0; i< node_path.size(); i++ )
        printNode ( node_path[i] );
      return  false;
    }
    reroute.erase ( edge_id );
    detour_cost -= E[edge_id].weight;
    i++;
  }
  i=node_path.size()-1;
  while ( i-1>=0 && isDanglingNode(bnet, node_path[i], reroute )) {
    int edge_id = getEdgeIdWithTwoEndingPoints ( node_path[i], node_path[i-1] );
    if ( edge_id==-1 ) {
      printNode(node_path[i]);
      printNode(node_path[i-1]);
      cout<< "3" << endl;
      for ( size_t i=0; i< node_path.size(); i++ )
        printNode ( node_path[i] );
      return  false;
    }
    reroute.erase ( edge_id );
    detour_cost -= E[edge_id].weight;
    //    printEdge ( edge_id );
    //    cout << "-" << E[edge_id].weight << endl;
    i--;
  }
  return  true;
  //  cout << "after 2 " << detour_cost << endl;
}

void  CGR::detourSingleEdge ( GridBox &box, int bnet_index, int edge_id )
{
  CoarseNet &bnet = bnets[bnet_index];
  int s = E[edge_id].srcPoint;
  int t = E[edge_id].snkPoint;
  assert ( E[edge_id].type!="VIA" );
  //  only consider 6 candidate detouring
  double  best_detour_cost = 100000;  //  infinity
  bnet.edges.erase ( edge_id );
  set<int>  best_reroute;
  //  horizontal edge
  if ( V[s].x!=V[t].x ) {
    assert ( s+1==t );
    // left
    if ( V[s].y+1 < fullGridY ) {
      if ( V[s].z+1 < numLayers ) {
        double  detour_cost = -E[edge_id].weight;
        double  &ref_detour_cost = detour_cost;
        set<int>  detour_route = bnet.edges;
        set<int>  &reroute = detour_route;
        int arr[] = {
          s, s+fullGridX*fullGridY, s+fullGridX*fullGridY+fullGridX, s+fullGridX,
          t+fullGridX, t+fullGridX*fullGridY+fullGridX, t+fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        if(!detourPath ( bnet, path, ref_detour_cost, reroute ) ) {
          printEdge ( edge_id );
          exit(0);
        }
        if ( detour_cost <=  best_detour_cost ) {
          best_detour_cost = detour_cost;
          best_reroute = detour_route;
        }
      }
      if ( V[s].z-1 > 0 ) {
        double  detour_cost = -E[edge_id].weight;
        double  &ref_detour_cost = detour_cost;
        set<int>  detour_route = bnet.edges;
        set<int>  &reroute = detour_route;
        int arr[] = {
          s, s-fullGridX*fullGridY, s-fullGridX*fullGridY+fullGridX, s+fullGridX,
          t+fullGridX, t-fullGridX*fullGridY+fullGridX, t-fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        if(!detourPath ( bnet, path, ref_detour_cost, reroute ) ) {
          printEdge ( edge_id );
          exit(0);
        }
        if ( detour_cost <=  best_detour_cost ) {
          best_detour_cost = detour_cost;
          best_reroute = detour_route;
        }
      }
    }
    //  right
    if ( V[s].y-1 >= 0 ) {
      if ( V[s].z+1 < numLayers ) {
        double  detour_cost = -E[edge_id].weight;
        double  &ref_detour_cost = detour_cost;
        set<int>  detour_route = bnet.edges;
        set<int>  &reroute = detour_route;
        int arr[] = {
          s, s+fullGridX*fullGridY, s+fullGridX*fullGridY-fullGridX, s-fullGridX,
          t-fullGridX, t+fullGridX*fullGridY-fullGridX, t+fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        if(!detourPath ( bnet, path, ref_detour_cost, reroute ) ) {
          printEdge ( edge_id );
          exit(0);
        }
        if ( detour_cost <=  best_detour_cost ) {
          best_detour_cost = detour_cost;
          best_reroute = detour_route;
        }
      }
      if ( V[s].z-1 > 0 ) {
        double  detour_cost = -E[edge_id].weight;
        double  &ref_detour_cost = detour_cost;
        set<int>  detour_route = bnet.edges;
        set<int>  &reroute = detour_route;
        int arr[] = {
          s, s-fullGridX*fullGridY, s-fullGridX*fullGridY-fullGridX, s-fullGridX,
          t-fullGridX, t-fullGridX*fullGridY-fullGridX, t-fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        if(!detourPath ( bnet, path, ref_detour_cost, reroute ) ) {
          printEdge ( edge_id );
          exit(0);
        }
        if ( detour_cost <=  best_detour_cost ) {
          best_detour_cost = detour_cost;
          best_reroute = detour_route;
        }
      }
    }
  }
  if ( V[s].y!=V[t].y ) {
    assert ( s+fullGridX==t );
    // left
    if ( V[s].x+1 < fullGridX ) {
      if ( V[s].z+1 < numLayers ) {
        double  detour_cost = -E[edge_id].weight;
        double  &ref_detour_cost = detour_cost;
        set<int>  detour_route = bnet.edges;
        set<int>  &reroute = detour_route;
        int arr[] = {
          s, s+fullGridX*fullGridY, s+fullGridX*fullGridY+1, s+1,
          t+1, t+fullGridX*fullGridY+1, t+fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        if(!detourPath ( bnet, path, ref_detour_cost, reroute ) ) {
          printEdge ( edge_id );
          exit(0);
        }
        if ( detour_cost <=  best_detour_cost ) {
          best_detour_cost = detour_cost;
          best_reroute = detour_route;
        }
      }
      if ( V[s].z-1 > 0 ) {
        double  detour_cost = -E[edge_id].weight;
        double  &ref_detour_cost = detour_cost;
        set<int>  detour_route = bnet.edges;
        set<int>  &reroute = detour_route;
        int arr[] = {
          s, s-fullGridX*fullGridY, s-fullGridX*fullGridY+1, s+1,
          t+1, t-fullGridX*fullGridY+1, t-fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        if(!detourPath ( bnet, path, ref_detour_cost, reroute ) ) {
          printEdge ( edge_id );
          exit(0);
        }
        if ( detour_cost <=  best_detour_cost ) {
          best_detour_cost = detour_cost;
          best_reroute = detour_route;
        }
      }
    }
    //  right
    if ( V[s].x-1 >= 0 ) {
      if ( V[s].z+1 < numLayers ) {
        double  detour_cost = -E[edge_id].weight;
        double  &ref_detour_cost = detour_cost;
        set<int>  detour_route = bnet.edges;
        set<int>  &reroute = detour_route;
        int arr[] = {
          s, s+fullGridX*fullGridY, s+fullGridX*fullGridY-1, s-1,
          t-1, t+fullGridX*fullGridY-1, t+fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        if(!detourPath ( bnet, path, ref_detour_cost, reroute ) ) {
          printEdge ( edge_id );
          exit(0);
        }
        if ( detour_cost <=  best_detour_cost ) {
          best_detour_cost = detour_cost;
          best_reroute = detour_route;
        }
      }
      if ( V[s].z-1 > 0 ) {
        double  detour_cost = -E[edge_id].weight;
        double  &ref_detour_cost = detour_cost;
        set<int>  detour_route = bnet.edges;
        set<int>  &reroute = detour_route;
        int arr[] = {
          s, s-fullGridX*fullGridY, s-fullGridX*fullGridY-1, s-1,
          t-1, t-fullGridX*fullGridY-1, t-fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        if(!detourPath ( bnet, path, ref_detour_cost, reroute ) ) {
          printEdge ( edge_id );
          exit(0);
        }
        if ( detour_cost <=  best_detour_cost ) {
          best_detour_cost = detour_cost;
          best_reroute = detour_route;
        }
      }
    }
  }

  //  vertical edge
  if ( V[s].z-2 > 0 ) {
    double  detour_cost = -E[edge_id].weight;
    double  &ref_detour_cost = detour_cost;
    set<int>  detour_route = bnet.edges;
    set<int>  &reroute = detour_route;
    int arr[] = {
      s, s-fullGridX*fullGridY, s-2*fullGridX*fullGridY, 
      t-2*fullGridX*fullGridY, t-fullGridX*fullGridY, t
    };
    vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
    if(!detourPath ( bnet, path, ref_detour_cost, reroute ) ) {
      printEdge ( edge_id );
      exit(0);
    }
    if ( detour_cost <=  best_detour_cost ) {
      best_detour_cost = detour_cost;
      best_reroute = detour_route;
    }
  }
  if ( V[s].z+2 < numLayers ) {
    double  detour_cost = -E[edge_id].weight;
    double  &ref_detour_cost = detour_cost;
    set<int>  detour_route = bnet.edges;
    set<int>  &reroute = detour_route;
    int arr[] = {
      s, s+fullGridX*fullGridY, s+2*fullGridX*fullGridY, 
      t+2*fullGridX*fullGridY, t+fullGridX*fullGridY, t
    };
    vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
    if(!detourPath ( bnet, path, ref_detour_cost, reroute ) ) {
      printEdge ( edge_id );
      exit(0);
    }
    if ( detour_cost <=  best_detour_cost ) {
      best_detour_cost = detour_cost;
      best_reroute = detour_route;
    }
  }
  if ( best_detour_cost <= 0 ) {
    bnet.edges = best_reroute;
  }
}

bool  CGR::isValidRoute ( set<int> &terminal_pins ) 
{
  //  dfs
  //  check if all the terminal pins are connected in cgr.graph
  unvisited_nodes = nodes;
  unvisited_terminal_pins = terminal_pins;
  dfs( *(unvisited_nodes.begin()) );
  return  unvisited_terminal_pins.empty() ;
}

void  CGR::dfs ( int u )
{
  unvisited_terminal_pins.erase(u);
  unvisited_nodes.erase(u);
  //  dfs recursively for each unvisited neighboring node
  for ( set<int>::iterator itr = graph[u].begin(); itr!= graph[u].end(); itr++ ) {
    //  unvisited
    if ( unvisited_nodes.find(*itr)==unvisited_nodes.end() )
      continue;
    //  recursively
    dfs ( *itr );
  }
}

//  add new candidate route
//    based on the route stored in bnet.edges
void  CGR::genCandidateRoute3 ( GridBox &box, int bnet_index, int M3 )
{
  bool  debug_mode = false;

  CoarseNet &bnet = bnets[bnet_index];
  bnet.congested_edges.clear();

  //  debug
  if ( debug_mode )
    for ( set<int>::iterator itr = bnet.edges.begin(); itr!= bnet.edges.end(); itr++ ) 
      assert ( E[*itr].cap!=0 );

  //  if the net is simple
  //    then use simple routing
  if ( bnet.terminal_pins.size()==2 && getNetHPWL(bnet)<20 ) {
    genCandidateRoutingForSimpleNet ( box, bnet, 2 );
    //    return;
  }

  //  identify all the congested edges in current route

  for ( set<int>::iterator itr = bnet.edges.begin(); itr != bnet.edges.end(); itr++ ) 
    if ( E[*itr].congestion_ratio > TCR ) 
      bnet.congested_edges.insert ( *itr );

  cout << "bnet " << bnet_index << " " << bnet.congested_edges.size() << " " << getNetHPWL ( bnet ) << endl;
  return;

  int original_wire_length = bnet.edges.size();
  //  create bundles
  for ( set<int>::iterator itr = bnet.congested_edges.begin(); itr!= bnet.congested_edges.end(); itr++ ) {
    if ( bnet.edges.size() > 2 * original_wire_length )
      break;
    createBundles ( box, bnet_index, *itr );
  }
  //  debug
  if ( debug_mode )
    for ( set<int>::iterator itr = bnet.edges.begin(); itr!= bnet.edges.end(); itr++ ) 
      assert ( E[*itr].cap!=0 );

  //  create pattern routing
  createPatternRoute ( box, bnet_index );

  buildGraph ( bnet.terminal_pins, bnet.edges );

  //  delete congested edge on this net
  for ( set<int>::iterator itr = bnet.congested_edges.begin(); itr!= bnet.congested_edges.end(); itr++ ) {
    bnet.edges.erase( *itr );
    deleteEdgeFromGraph ( bnet.terminal_pins, *itr );
    if ( isValidRoute ( bnet.terminal_pins )==false ) {
      bnet.edges.insert ( *itr );
      addEdgeToGraph ( *itr );
    }
  }

  //  sort edge weights
  vector<pair<int, double> > edge_weight;
  for ( set<int>::iterator itr = bnet.edges.begin(); itr!= bnet.edges.end(); itr++ ) {
    assert ( E[*itr].cap!=0 );
    edge_weight.push_back ( make_pair(*itr, E[*itr].weight) );
  }
  sort ( edge_weight.begin(), edge_weight.end(), compare_pair_int_double );

  //  build graph
  buildGraph ( bnet.terminal_pins, bnet.edges );
  assert ( isValidRoute ( bnet.terminal_pins ) );

  //  reverse-delete
  for ( vector<pair<int,double> >::iterator itr = edge_weight.begin(); itr!=edge_weight.end(); itr++ ) {
    bnet.edges.erase( (*itr).first );
    deleteEdgeFromGraph ( bnet.terminal_pins, (*itr).first );
    if ( isValidRoute ( bnet.terminal_pins )==false ) {
      bnet.edges.insert ( (*itr).first );
      addEdgeToGraph ( (*itr).first );
    }
  }

  for ( set<int>::iterator it = bnet.edges.begin(); it!= bnet.edges.end(); it++ ) {
    //    E[*it].weight += DELTA_EDGE_WEIGHT;
    if ( E[*it].cap==0 || map_from_edge_id_to_index_in_E_id.count(*it) == 0 ) {
      cout << "ERROR" << endl;
      printEdge(*it);
      exit(0);
    }
  }

  NetRoute new_candidate_route;
  new_candidate_route.edges = bnet.edges;
  bnet.candidate_routes.push_back ( new_candidate_route );

  //
  //  print result
  if ( debug_mode ) {
    string  str = bnet.candidate_routes[0].edges!=bnet.candidate_routes[1].edges?"TRUE":"FALSE" ;
    cout << bnet.old_index << " " << bnet.terminal_pins.size() << " " << getNetHPWL ( bnet ) << " " << bnet.congested_edges.size() << " " << bnet.candidate_routes.size() << " " << str << " " << bnet.candidate_routes[0].edges.size() << " " << bnet.candidate_routes[1].edges.size() << endl;
    cout << "route 0" << endl;
    for ( set<int>::iterator it = bnet.candidate_routes[0].edges.begin(); it!=bnet.candidate_routes[0].edges.end(); it++ )
      printEdge(*it);
    cout << "route 1" << endl;
    for ( set<int>::iterator it = bnet.candidate_routes[1].edges.begin(); it!=bnet.candidate_routes[1].edges.end(); it++ )
      printEdge(*it);
  }
}

void  CGR::createBundles ( GridBox &box, int bnet_index, int edge_id )
{
  //  cout << endl<< " # congested edge : " ;
  //  printEdge ( edge_id );
  CoarseNet &bnet = bnets[bnet_index];
  int s = E[edge_id].srcPoint;
  int t = E[edge_id].snkPoint;
  if ( E[edge_id].type=="VIA" ) {
    return;
  }
  //  horizontal edge
  if ( V[s].x!=V[t].x ) {
    assert ( s+1==t );
    // left
    if ( V[s].y+1<box.y0+box.dy && E[edge_id+fullGridX-1].congestion_ratio <= E[edge_id].congestion_ratio ) {
      if ( V[s].z+1 < numLayers ) {
        int arr[] = {
          s, s+fullGridX*fullGridY, s+fullGridX*fullGridY+fullGridX, s+fullGridX,
          t+fullGridX, t+fullGridX*fullGridY+fullGridX, t+fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        for ( size_t i=0; i < path.size()-1; i++ ) {
          int edge_id = getEdgeIdWithTwoEndingPoints ( path[i], path[i+1] );
          assert ( E[edge_id].cap!=0 );
          bnet.edges.insert ( edge_id );
        }
      }
      if ( V[s].z-1 > 0 ) {
        int arr[] = {
          s, s-fullGridX*fullGridY, s-fullGridX*fullGridY+fullGridX, s+fullGridX,
          t+fullGridX, t-fullGridX*fullGridY+fullGridX, t-fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        for ( size_t i=0; i < path.size()-1; i++ ) {
          int edge_id = getEdgeIdWithTwoEndingPoints ( path[i], path[i+1] );
          assert ( E[edge_id].cap!=0 );
          bnet.edges.insert ( edge_id );
        }
      }
    }
    //  right
    if ( V[s].y-1>=box.y0 && E[edge_id-fullGridX+1].congestion_ratio <= E[edge_id].congestion_ratio ) {
      if ( V[s].z+1 < numLayers ) {
        int arr[] = {
          s, s+fullGridX*fullGridY, s+fullGridX*fullGridY-fullGridX, s-fullGridX,
          t-fullGridX, t+fullGridX*fullGridY-fullGridX, t+fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        for ( size_t i=0; i < path.size()-1; i++ ) {
          int edge_id = getEdgeIdWithTwoEndingPoints ( path[i], path[i+1] );
          assert ( E[edge_id].cap!=0 );
          bnet.edges.insert ( edge_id );
        }
      }
      if ( V[s].z-1 > 0 ) {
        int arr[] = {
          s, s-fullGridX*fullGridY, s-fullGridX*fullGridY-fullGridX, s-fullGridX,
          t-fullGridX, t-fullGridX*fullGridY-fullGridX, t-fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        for ( size_t i=0; i < path.size()-1; i++ ) {
          int edge_id = getEdgeIdWithTwoEndingPoints ( path[i], path[i+1] );
          assert ( E[edge_id].cap!=0 );
          bnet.edges.insert ( edge_id );
        }
      }
    }
  }
  if ( V[s].y!=V[t].y ) {
    assert ( s+fullGridX==t );
    // left
    if ( V[s].x+1<box.x0+box.dx && E[edge_id+fullGridY-1].congestion_ratio <= E[edge_id].congestion_ratio ) {
      if ( V[s].z+1 < numLayers ) {
        int arr[] = {
          s, s+fullGridX*fullGridY, s+fullGridX*fullGridY+1, s+1,
          t+1, t+fullGridX*fullGridY+1, t+fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        for ( size_t i=0; i < path.size()-1; i++ ) {
          int edge_id = getEdgeIdWithTwoEndingPoints ( path[i], path[i+1] );
          assert ( E[edge_id].cap!=0 );
          bnet.edges.insert ( edge_id );
        }
      }
      if ( V[s].z-1 > 0 ) {
        int arr[] = {
          s, s-fullGridX*fullGridY, s-fullGridX*fullGridY+1, s+1,
          t+1, t-fullGridX*fullGridY+1, t-fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        for ( size_t i=0; i < path.size()-1; i++ ) {
          int edge_id = getEdgeIdWithTwoEndingPoints ( path[i], path[i+1] );
          assert ( E[edge_id].cap!=0 );
          bnet.edges.insert ( edge_id );
        }
      }
    }
    //  right
    if ( V[s].x-1>=box.x0 && E[edge_id-fullGridY+1].congestion_ratio <= E[edge_id].congestion_ratio ) {
      if ( V[s].z+1 < numLayers ) {
        int arr[] = {
          s, s+fullGridX*fullGridY, s+fullGridX*fullGridY-1, s-1,
          t-1, t+fullGridX*fullGridY-1, t+fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        for ( size_t i=0; i < path.size()-1; i++ ) {
          int edge_id = getEdgeIdWithTwoEndingPoints ( path[i], path[i+1] );
          assert ( E[edge_id].cap!=0 );
          bnet.edges.insert ( edge_id );
        }
      }
      if ( V[s].z-1 > 0 ) {
        int arr[] = {
          s, s-fullGridX*fullGridY, s-fullGridX*fullGridY-1, s-1,
          t-1, t-fullGridX*fullGridY-1, t-fullGridX*fullGridY, t
        };
        vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
        for ( size_t i=0; i < path.size()-1; i++ ) {
          int edge_id = getEdgeIdWithTwoEndingPoints ( path[i], path[i+1] );
          assert ( E[edge_id].cap!=0 );
          bnet.edges.insert ( edge_id );
        }
      }
    }
  }

  //  vertical edge
  if ( V[s].z-2>0 && E[edge_id - 2*edgeCount2D].congestion_ratio <= E[edge_id].congestion_ratio )  {
    int arr[] = {
      s, s-fullGridX*fullGridY, s-2*fullGridX*fullGridY, 
      t-2*fullGridX*fullGridY, t-fullGridX*fullGridY, t
    };
    vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );
    for ( size_t i=0; i < path.size()-1; i++ ) {
      int edge_id = getEdgeIdWithTwoEndingPoints ( path[i], path[i+1] );
      bnet.edges.insert ( edge_id );
    }
  }
  if ( V[s].z+2<numLayers && E[edge_id + 2*edgeCount2D].congestion_ratio <= E[edge_id].congestion_ratio )  {
    int arr[] = {
      s, s+fullGridX*fullGridY, s+2*fullGridX*fullGridY, 
      t+2*fullGridX*fullGridY, t+fullGridX*fullGridY, t
    };
    vector<int> path ( arr, arr+sizeof(arr)/sizeof(int) );

    for ( size_t i=0; i < path.size()-1; i++ ) {
      int edge_id = getEdgeIdWithTwoEndingPoints ( path[i], path[i+1] );
      bnet.edges.insert ( edge_id );
    }
  }
}

void  CGR::deleteEdgeFromGraph ( set<int> &terminal_pins, int edge_id )
{
  int u = E[edge_id].srcPoint;
  int v = E[edge_id].snkPoint;
  assert ( graph[u].find(v)!=graph[u].end() );
  graph[u].erase(v);
  //  if ( graph[u].empty() && terminal_pins.find(u)==terminal_pins.end() ) {
  //    nodes.erase(u);
  //  }
  assert ( graph[v].find(u)!=graph[v].end() );
  graph[v].erase(u);
  //  if ( graph[v].empty() && terminal_pins.find(v)==terminal_pins.end() ) {
  //    nodes.erase(v);
  //  }
}

void  CGR::addEdgeToGraph ( int edge_id ) 
{
  int u = E[edge_id].srcPoint;
  int v = E[edge_id].snkPoint;
  nodes.insert ( u );
  nodes.insert ( v );
  assert ( graph.count(u)!=0 );
  graph[u].insert (v);
  assert ( graph.count(v)!=0 );
  graph[v].insert (u);
}

void  CGR::buildGraph ( set<int> &pins, set<int> &edges )
{
  graph.clear();
  nodes.clear();
  for ( set<int>::iterator itr = edges.begin(); itr != edges.end(); itr++ ) {
    int u = E[*itr].srcPoint;
    int v = E[*itr].snkPoint;
    nodes.insert ( u );
    nodes.insert ( v );

    //
    if ( graph.count(u)==0 ) {
      set<int>  tmp;
      tmp.insert ( v );
      graph.insert ( make_pair(u,tmp) );
    }
    else
      graph[u].insert (v);

    //
    if ( graph.count(v)==0 ) {
      set<int>  tmp;
      tmp.insert ( u );
      graph.insert ( make_pair(v,tmp) );
    }
    else
      graph[v].insert (u);

  }
}

void  CGR::createPatternRoute ( GridBox &box, int bnet_index )
{
  CoarseNet &bnet = bnets[bnet_index];
  for ( int i=0; i < bnet.terminal_pins.size()-1; i++ ) {
    for ( set<int>::iterator it = bnet.subnets[i].edges.begin(); it!=bnet.subnets[i].edges.end(); it++ )
      if ( E[*it].congestion_ratio > TCR ) {
        addSimplePattern ( bnet.subnets[i].s, bnet.subnets[i].t, bnet_index );
        break;
      }
  }
}

void  CGR::addSimplePattern ( int s, int t, int bnet_index )
{
  set<int> &edges = bnets[bnet_index].edges;
  if ( V[s].z==V[t].z ) {
    if ( V[s].x!=V[t].x && V[s].y==V[t].y  ) {
      if ( V[s].z+2<numLayers && horizCaps[V[s].z+2]!=0 ) {
        set<int>  seg1 = getFlatSegment ( s, s+2*fullGridX*fullGridY );
        set<int>  seg2 = getFlatSegment ( s+2*fullGridX*fullGridY, t+2*fullGridX*fullGridY );
        set<int>  seg3 = getFlatSegment ( t, t+2*fullGridX*fullGridY ); 
        edges.insert ( seg1.begin(), seg1.end() );
        edges.insert ( seg2.begin(), seg2.end() );
        edges.insert ( seg3.begin(), seg3.end() );
      }
      else {
        if ( V[s].z+1<numLayers && horizCaps[V[s].z+1]!=0 ) {
          set<int>  seg1 = getFlatSegment ( s, s+fullGridX*fullGridY );
          set<int>  seg2 = getFlatSegment ( s+fullGridX*fullGridY, t+fullGridX*fullGridY );
          set<int>  seg3 = getFlatSegment ( t, t+fullGridX*fullGridY ); 
          edges.insert ( seg1.begin(), seg1.end() );
          edges.insert ( seg2.begin(), seg2.end() );
          edges.insert ( seg3.begin(), seg3.end() );
        }
      }
    }
    else  if ( V[s].y!=V[t].y && V[s].x==V[t].x ) {
      if ( V[s].z+2<numLayers && vertCaps[V[s].z+2]!=0 ) {
        set<int>  seg1 = getFlatSegment ( s, s+2*fullGridX*fullGridY );
        set<int>  seg2 = getFlatSegment ( s+2*fullGridX*fullGridY, t+2*fullGridX*fullGridY );
        set<int>  seg3 = getFlatSegment ( t, t+2*fullGridX*fullGridY ); 
        edges.insert ( seg1.begin(), seg1.end() );
        edges.insert ( seg2.begin(), seg2.end() );
        edges.insert ( seg3.begin(), seg3.end() );
      }
      else {
        if ( V[s].z+1<numLayers && vertCaps[V[s].z+1]!=0 ) {
          set<int>  seg1 = getFlatSegment ( s, s+fullGridX*fullGridY );
          set<int>  seg2 = getFlatSegment ( s+fullGridX*fullGridY, t+fullGridX*fullGridY );
          set<int>  seg3 = getFlatSegment ( t, t+fullGridX*fullGridY ); 
          edges.insert ( seg1.begin(), seg1.end() );
          edges.insert ( seg2.begin(), seg2.end() );
          edges.insert ( seg3.begin(), seg3.end() );
        }
      }
    }
    else {
      //  L shape
      int u = xyzToGridNetTerminalId ( V[s].x, V[t].y, V[s].z );
      int v = xyzToGridNetTerminalId ( V[t].x, V[s].y, V[s].z );
      addSimplePattern ( s, u, bnet_index );
      addSimplePattern ( u, t, bnet_index );
      addSimplePattern ( s, v, bnet_index );
      addSimplePattern ( v, t, bnet_index );
    }
  }
}

void  CGR::getInputGlobalRoutingInfoInBox ( std::ofstream &stream_file ) 
{
  stream_file << "local metrics : ";
  stream_file << getACEInBox ( 0.001 ) << " ";
  stream_file << getACEInBox ( 0.005 ) << " ";
  stream_file << getACEInBox ( 0.01 ) << " ";
  stream_file << getACEInBox ( 0.02 ) << " ";
  stream_file << getACEInBox ( 0.05 ) << " ";
  stream_file << getACEInBox ( 0.1 ) << " ";
  stream_file << getACEInBox ( 0.2 ) << " ";
  stream_file << endl;
}

double  CGR::getACEInBox ( double y ) 
{
  vector<double>  edge_congestion_ratio;
  for ( vector<int>::iterator itr_edge = E_id.begin(); itr_edge != E_id.end(); itr_edge++ ) 
    if ( E[*itr_edge].cap!=0 && E[*itr_edge].type!="VIA" )
      edge_congestion_ratio.push_back ( E[*itr_edge].congestion_ratio );
  sort ( edge_congestion_ratio.begin(), edge_congestion_ratio.end() );
  reverse ( edge_congestion_ratio.begin(), edge_congestion_ratio.end() );

  double  sum_congestion_edge = 0;
  int edge_count = 0;
  for ( size_t i=0; i < edge_congestion_ratio.size() * y ; i++ ) {
    sum_congestion_edge += edge_congestion_ratio[i];
    ++ edge_count;
  }
//  cout << "SMALLEST ("<<y<<") = " << edge_congestion_ratio[edge_congestion_ratio.size()*y] << endl;
  return  (double)sum_congestion_edge / edge_count;
}

int CGR::getTWL (void)
{
  int cnt = 0;
  for ( size_t i=0; i < nets.size(); i++ )
    cnt+= nets[i].edges.size();
  return  cnt;
}

double  CGR::getOBJInBox (void) 
{
  double  obj = 0;
  for ( vector<int>::iterator it = E_id.begin(); it != E_id.end(); it++ )
    if ( E[*it].cap!=0 && E[*it].type!="VIA" )
      obj += getEdgePenalty(*it);
  return  obj;
}

void  CGR::separateObjVal (void)
{
  obj_good = 0;
  obj_bad = 0;
  for ( vector<int>::iterator itr = E_id.begin(); itr!= E_id.end(); itr++ ) {
    if ( E[*itr].cap==0 || E[*itr].type=="VIA" )
      continue;
    if ( E[*itr].congestion_ratio > TCR )
      obj_bad += getEdgePenalty ( *itr );
    else
      obj_good += getEdgePenalty ( *itr );
  }
  cout << "OBJVAL_GOOD = " << obj_good << endl;
  cout << "OBJVAL_BAD = " << obj_bad << endl;
  cout << "OBJVAL = " << obj_good+obj_bad << endl;
}

void  CGR::computeBoxObjVal (void)
{
  obj_bad = 0;
  for ( vector<int>::iterator itr = E_id.begin(); itr!= E_id.end(); itr++ ) {
    if ( E[*itr].cap==0 || E[*itr].type=="VIA" )
      continue;
    if ( E[*itr].congestion_ratio>TCR )
      obj_bad += getEdgePenalty ( *itr );
//      obj_bad += slope * (E[*itr].congestion_ratio-TCR);
  }
}

void  CGR::recoverNets (void)
{
  for ( size_t i=0; i < bnets.size(); i++ )
  {
    CoarseNet &bnet = bnets[i];
    //  rip-up
    for ( set<int>::iterator itr = bnet.edges.begin(); itr!= bnet.edges.end(); itr++ ) {
      assert ( find(E_id.begin(), E_id.end(), *itr )!=E_id.end() );
      assert ( nets[bnet.old_index].edges.find(*itr)!=nets[bnet.old_index].edges.end() );
      nets[ bnet.old_index ].edges.erase ( *itr );
      assert ( E[*itr].nets_passing_through.find(bnet.old_index) != E[*itr].nets_passing_through.end() );
      E[*itr].nets_passing_through.erase ( bnet.old_index );
      -- E[*itr].util;
    }
    //  check if the remained edges are not in box
    for ( set<int>::iterator itr = nets[bnet.old_index].edges.begin(); itr!= nets[bnet.old_index].edges.end(); itr++ ) {
      assert ( find(E_id.begin(), E_id.end(), *itr )==E_id.end() );
    }
    //  recover
    bnet.edges = bnet.edges_backup;
    //  reroute
    for ( set<int>::iterator itr = bnet.edges.begin(); itr!= bnet.edges.end(); itr++ ) {
      assert ( find(E_id.begin(), E_id.end(), *itr )!=E_id.end() );
      assert ( nets[bnet.old_index].edges.find(*itr)==nets[bnet.old_index].edges.end() );
      nets[ bnet.old_index ].edges.insert ( *itr );
      assert ( E[*itr].nets_passing_through.find(bnet.old_index) == E[*itr].nets_passing_through.end() );
      E[*itr].nets_passing_through.insert ( bnet.old_index );
      ++ E[*itr].util;
    }
    generateLongestSegments ( bnet.old_index );
  }
  for ( vector<int>::iterator it = E_id.begin(); it!= E_id.end(); it++ )
    E[*it].congestion_ratio = (double)(E[*it].blk + E[*it].util)/(E[*it].cap+0.0001);
}

void  CGR::genCandidateRoute4 ( GridBox &box, int bnet_index, int M3 )
{
  bool  debug_mode = false;

  CoarseNet &bnet = bnets[bnet_index];

  bnet.congested_edges.clear();

  //  if the net is simple
  //    then use simple routing
  if ( bnet.terminal_pins.size()==2 && getNetHPWL(bnet)<20 ) {
    genCandidateRoutingForSimpleNet ( box, bnet, M3 );
    return;
  }

  //  identify all the congested edges in current route
  for ( set<int>::iterator itr = bnet.edges.begin(); itr != bnet.edges.end(); itr++ ) 
    if ( E[*itr].congestion_ratio > TCR ) 
      bnet.congested_edges.insert ( *itr );

  int original_wire_length = bnet.edges.size();
  //  create bundles
  for ( set<int>::iterator itr = bnet.congested_edges.begin(); itr!= bnet.congested_edges.end(); itr++ ) {
    if ( bnet.edges.size() > 2 * original_wire_length )
      break;
    createBundles ( box, bnet_index, *itr );
  }

  //  create pattern routing
  createPatternRoute ( box, bnet_index );

  //  sort edge weights
  vector<pair<int, double> > edge_weight;
  for ( set<int>::iterator itr = bnet.edges.begin(); itr!= bnet.edges.end(); itr++ ) {
    assert ( E[*itr].cap!=0 );
    edge_weight.push_back ( make_pair(*itr, E[*itr].weight) );
  }
  sort ( edge_weight.begin(), edge_weight.end(), compare_pair_int_double );


  set<int>  edges_backup = bnet.edges;
  //  generate M3 candidate routes using reverse-delete
  for ( int ii =0, jj=0; ii< 2*M3 && ii<edges_backup.size() && jj< M3; ii++ ) {
    //  build graph
    bnet.edges = edges_backup;
    buildGraph ( bnet.terminal_pins, edges_backup );
    assert ( isValidRoute ( bnet.terminal_pins ) );
    //  reverse-delete
    for ( vector<pair<int,double> >::iterator itr = edge_weight.begin() + ii; itr!=edge_weight.end(); itr++ ) {
      bnet.edges.erase( (*itr).first );
      deleteEdgeFromGraph ( bnet.terminal_pins, (*itr).first );
      if ( isValidRoute ( bnet.terminal_pins )==false ) {
        bnet.edges.insert ( (*itr).first );
        addEdgeToGraph ( (*itr).first );
      }
    }
    //  check if this is already in Ti
    bool  isIncluded = false;
    for ( size_t j=0; j < bnet.candidate_routes.size(); j++ )
      if ( bnet.candidate_routes[j].edges == bnet.edges ) {
        isIncluded = true;
        break;
      }
    if ( isIncluded == false ) {
      NetRoute new_candidate_route;
      new_candidate_route.edges = bnet.edges;
      bnet.candidate_routes.push_back ( new_candidate_route );
      ++ jj;
    }
  }
}

vector<double>  CGR::getHist ( std::ofstream &stream_file )
{
  int hist_num = 31;
  bool  debug_mode = true;
  vector<double>  hist(hist_num, 0);
  int cnt = 0;
  for ( vector<int>::iterator it= E_id.begin(); it!= E_id.end(); it++ ) {
    if ( E[*it].cap==0 || E[*it].type=="VIA" )
      continue;
    int key = E[*it].congestion_ratio / 0.1;
    ++ cnt;
    if ( key>=hist_num )
      key = hist_num-1;
    ++ hist[key];
  }
  for ( size_t i=0; i < hist_num; i++ )
    hist[i] = hist[i] / cnt;

  if ( debug_mode ) {
    stream_file << "congestion ratio distribution : " << endl;
    for ( size_t i=0; i<hist_num; i++ )
      stream_file << hist[i] << " ";
    stream_file << endl;
  }

  return  hist;
}

double  CGR::getObjValInBox ( void )
{
  double  objval =0;
  //  part 1
  for ( vector<int>::iterator it = E_id.begin(); it!= E_id.end(); it++ ) {
    E[*it].virtual_util = E[*it].util;
  }
  //  part 2
  for ( vector<CoarseNet>::iterator it = bnets.begin(); it!= bnets.end(); it++ )
    for ( set<int>::iterator itr_edge = (*it).edges.begin(); itr_edge!= (*it).edges.end(); itr_edge++ ) {
      assert ( E[*itr_edge].cap!=0 );
      ++ E[*itr_edge].virtual_util;
    }
  //  part 3
  for ( vector<int>::iterator it = E_id.begin(); it!= E_id.end(); it++ ) {
    if ( (double)(E[*it].blk+E[*it].virtual_util)/(E[*it].cap+0.0001) > TCR )
//      objval += slope * ((double)(E[*it].blk+E[*it].virtual_util)/(E[*it].cap+0.0001) - TCR ); 
      objval+= getVirtualEdgePenalty ( *it );
  }
  return  objval;
}

double  CGR::getTHCR (double yy) 
{
  vector<double>  crv;
  for ( vector<int>::iterator itr = E_id.begin(); itr!= E_id.end(); itr++ )
    if ( E[*itr].cap != 0 && E[*itr].type!="VIA"  )
      crv.push_back ( E[*itr].congestion_ratio );
  sort ( crv.begin(), crv.end() );
  reverse ( crv.begin(), crv.end() );

  return crv[crv.size()*yy];
}

double  CGR::getVirtualEdgePenalty ( int edge_id )
{
  double  edge_penalty = 0;
  assert ( mj.size()==bj.size() );
  for ( size_t j=0; j < mj.size(); j++ ) {  //  get the value of edge penalty
    double  penalty_for_single_segment = mj[j] * ( E[edge_id].virtual_util +E[edge_id].blk ) / ( E[edge_id].cap+0.0001 ) + bj[j];
    edge_penalty = ( penalty_for_single_segment > edge_penalty ) ? penalty_for_single_segment : edge_penalty;
  } 
  return  edge_penalty;
}

double  CGR::getPenalty ( double y )
{
  double  edge_penalty = 0;
  assert ( mj.size()==bj.size() );
  for ( size_t j=0; j < mj.size(); j++ ) {  //  get the value of edge penalty
    double  penalty_for_single_segment = mj[j] * y  + bj[j];
    edge_penalty = ( penalty_for_single_segment > edge_penalty ) ? penalty_for_single_segment : edge_penalty;
  } 
  return  edge_penalty;
}

